{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the basic libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "import os\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the environment for running the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths for loading/storing results.\n",
    "name = 'results/Gaussian_MINEE' # filename\n",
    "chkpt_name = name+'.pt'              # checkpoint\n",
    "fig_name = name+'.pdf'               # output figure\n",
    "\n",
    "# use GPU if available\n",
    "if torch.cuda.is_available(): \n",
    "    torch.set_default_tensor_type(torch.cuda.FloatTensor)\n",
    "else:\n",
    "    torch.set_default_tensor_type(torch.FloatTensor)\n",
    "\n",
    "# initialize random seed\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data using the gaussian model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from data.gaussian import Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 400   # sample size\n",
    "rho = 0.9           # model parameter\n",
    "\n",
    "rep = 1             # number of repeated runs\n",
    "d = 6               # number of dimensions for X (and Y)\n",
    "\n",
    "X = np.zeros((rep,sample_size,d))\n",
    "Y = np.zeros((rep,sample_size,d))\n",
    "g = Gaussian(sample_size=sample_size,rho=rho)\n",
    "for i in range(rep):\n",
    "    for j in range(d):\n",
    "        data = g.data\n",
    "        X[i,:,j] = data[:,0]\n",
    "        Y[i,:,j] = data[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the first dimension of $Y$ against that of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XvQXVd53/HvIwsL40tAsRsH7FcmhpBQj4bLW4MnF2hwsGnd0CglJWEwmRQ05DYhQ20KJkC4ZIzd8ZAplKIWCgbCpcEEikkENE08NDaDxFBFYGhtwmsb2xgQjrBxbYSf/nH2kte7tO9nn7P2Oef3mdGg85599l7nFV7Puj7L3B0REZEtuQsgIiLjoIAgIiKAAoKIiBQUEEREBFBAEBGRggKCiIgACgiSgZn9tZm9eE7P+i0z+6aZ3WNmP9ri+q+b2fnzKNsYmNm7zeyNucsh46CAIDNRVKz3FRXxN83sv5rZSR3vcZaZuZlt7VmGhwFXAc9295Pc/Tt97lNzfzezxw15T5GcFBBklv6Fu58EPAX4J8Cr5/z8HwMeDnxpzs8VWUgKCDJz7v4N4C+Ac9L3zGyLmb3azDbM7C4zu9rMfqR4+7rif+8uehrnlXx+m5m9xcxuL/68pfjZTwJfjT7/V2VlM7MXFs/+jpldlrx3rpldb2Z3m9kdZvZWMzu+eC+U7X8XZfvXZvYoM/uEmX3LzL5b/P2Mqt+Lmb3CzL5hZt8zs6+a2bOanlu872b222b2f4vPvsHMzi4+c9jMPhyV85lmdpuZvcrMvl303F5QU6aLzOyLxbP/1sx2NpVXloi764/+DP4H+DpwfvH3M5m00t9QvP5r4MXF338TuAn4CeAk4BrgvcV7ZwEObK15zuuBG4B/BJwG/G30nNrPA08E7gF+HtjGZHjpSFTupwJPB7YW97oReFn0eQceF73+UeBXgEcAJwP/Dfjzimc/AbgVeHRU1rM7PPfjwCnAPwbuB/5H8Tv8EeDLwIuKa59ZfKeriu/4DOBe4AnF++8G3lj8/SnAXcDTgOOAFxX/jtvqyqs/y/NHPQSZpT83s7uBzwJ/A/xxyTUvAK5y96+5+z3AK4Hnd5g3eAHwene/y92/BfwR8MKWn/1XwCfc/Tp3vx/4Q+DB8Ka773f3G9z9iLt/HXgHkwq1lLt/x90/4u7fd/fvAW+quf6HTCraJ5rZw9z96+5+c4fnvtndD7v7l4CDwKeK3+E/MOmNPTm5/g/d/X53/xvgWuBXS8r0EuAd7v45d/+hu7+HSbB5el15ZXkoIMgs/Ut3f6S773D333b3+0queTSwEb3eYNIy/rGWzyj7/KM7fPbW8MLd7wWOTjyb2U8Wwz53mtlhJgHt1KqbmdkjzOwdxRDUYSZDXo80s+PSa939JuBlwOuAu8zsg2b26A7P/Wb09/tKXscT+N8tvltQ9TvaAby8GC66uwjmZzLpFVSWV5aHAoLkdjuTiihYYzLE8U0mQyN9Pn97y2ffwaTCAyYVOpNhn+DtwFeAx7v7KcCrAKu538uZDK08rbj+58Otyy529z91958tyu/Am3s+t8mjzOzE6HXV7+hW4E1FEA9/HuHuH2gorywJBQTJ7QPAH5jZY4tlqX8MfMjdjwDfYjKE8xMNn3+1mZ1mZqcCrwHe1/LZfwZcZGY/W0zCvp7N/02cDBwG7jGznwJ+K/n8N5OyncykdX63mW0HXlv1YDN7gpn9gpltA/5f8bkftnxuH39kZseb2c8BFzGZ30j9Z+ClZvY0mzjRzP65mZ3cUF5ZEgoIktu7gPcyGV75eyaVze8BuPv3mYzD/69iCOPpJZ9/I7APOAD8HfCF4meNivH33wH+lElv4bvAbdEl/xb4deB7TCrLDyW3eB3wnqJsvwq8BTgB+DaTie6/rHn8NuDy4to7mUyKv6rlc7u6k8l3ux14P/BSd/9KepG772Myj/DW4vqbgN9oUV5ZEuauA3JElpWZPRN4n7tXLn8VCdRDEBERQAFBREQK2YaMzOzhTMaNtzFZZvhn7l45CSciIrOVMyAYcKK732OTJGSfBX7f3W/IUiARkRXXK4vkEHwSie4pXj6s+FMbnU499VQ/66yzZlwyEZHlsn///m+7+2lN12ULCADFDs79wOOAt7n750qu2Q3sBlhbW2Pfvn3zLaSIyIIzs43mqzJPKhf5Up4EnAGca2bHZMN09z3uvu7u66ed1hjgRESkp1GsMnL3u5lkwLwwc1FERFZWtoBQpBp4ZPH3E4DzmeRvERGRDHLOIfw4k23/xzEJTB92909kLI+IyErLucroAMfmbBcRkUxGMYcgIiL5KSCIiIzYJVdfzyVXXz+XZykgiIgIkHljmoiIlAu9ggMbhza9vvLi82b2TPUQREQym+ewUB31EERERij0BObRMwgUEEREMskxLFRHAUFEZMTmGRwUEEREMskxLFRHk8oiIgMaywRxH+ohiIhklrtnECggiMhCG8twy9gmiPvQkJGIiABgk6ONF8P6+rrrCE0RgWNb5Dt3bAfyt8jH2DMws/3uvt50nXoIIiICqIcgIgtujC3ytuZVdvUQREQyW7QlqFplJCILbSw9gy6t/bGuSFJAEBEZ2Fgr/CYKCCIiUyir/G++8zBnn35K5WfGlrIiUEAQERnY2aefwpUXnze6Cr+JAoKIyBSmae2PLVAoIIiIzMjYKvwmCggiksWiDac0WYbvoX0IIiI9LNoegzbUQxCRuVqEJZljLNM8ZAsIZnYmcDVwOvAgsMfd/yRXeURE2liEgNZXzh7CEeDl7v4FMzsZ2G9mn3b3L2csk4jM2FjX4MO4Kvscz84WENz9DuCO4u/fM7MbgccACggiksXNdx5uvGbMAW1ao5hDMLOzgCcDnyt5bzewG2BtbW2u5RKR2RljRZruLs7ZM8jRS8keEMzsJOAjwMvc/Zjw7O57gD0wSX895+KJyApIK+ETtzVXjWMMaNPKGhDM7GFMgsH73f2anGURkeXXtrVdl4do2ns3yTkklXOVkQHvBG5096tylUNEZJnnBbrI2UP4GeCFwN+Z2ReLn73K3T+ZsUwisgSqNozNYlx+VmP+OYJSzlVGnwUs1/NFRFKr2jMIsk8qi4gMZdcVewG49/4jm36+c8f2Tf87i4p/lveeF+UyEhERQD0EEVkAZePy8c/C30PPICwbDauF+rba6+YD0rmDnTu2c/Odh7nk6usXtpegHoKItLaMGT7lIeohiMho1Z1XXLaqZ6gVPm1WDsVLVUPKi3vvP8KBjUMLu3xVAUFEGs1qaWXX+7TJNST9KSCIrLCxt2TTVng4vD78LL4mvT59r+9zm+5TNa+xiBQQRKTR0Dt52/Q44rmKm+88fMxwTJNFn+DNQQFBZAWNKe9/G/GcQZ30e00bFLp+bqy/v7YUEESktaEqvLoeR9lyzvh/08+m4nmGe+8/op5CBwoIIiuoyxBQ297DvHoZdb2bsNIoDDFBt8ylY+8pzZoCgohkU1bxNgWrpjmEOCjEk9DSzNwX58yZ9fV137dvX+5iiKyEqqGbsgo6bpGHXcLXXHrBIM+vqtCnfb/s2qbv2vW+Y2Fm+919vek69RBEZJRmnT56ESv2WVMPQURqNS0JhYda1lsMHiyqlLFn/2zKj1R1fZtexNiohyAiMxVW88ziuMlZtt4XbcntPCkgiCyooSuyqpxAdUtCw3xB3FqeZwU79LPa7Ehe5gCigCAinZSt84d2PYWm1nnX1nuf3EZ9KvZlDgIxBQSRBTP0kEd6v11X7D26Yqjt8s34/Xn2DOJyzuvZyxwUFBBEVsgQFWecRuLEbVs7rfVvap2XvR/OYIivTXsG0/QU6qzafIMCgsiC6Trk0XRdVSVcl+M/98avMDwVKuo+E9tyLAUEkRUQ8vmUHTTTt0Kf5rNNn4t7BnUpKrYYnHD81pntOViFieSYAoLIgmpbCcaZP9veL14t1OfZXdRVtk1lPuH4reodDEgBQWSBDNVSLTtfYF6VfFtVQ1NpoAs/a7OZrG8yv2XvGQQKCCIjNk3FWnbaWJszBco+37cMVcrKtOuKvUcr/qHPNZB2FBBEBjLLceZdV+wF2i2zLHsvrWi7biBrWyF3aZHf98CR1quD4mGhPsdbdukZrMqKojJZA4KZvQu4CLjL3c/JWRaRMUkrpy3W/15DJIlre2xlk/R73ffAkaP5j+JgN23vZBUr8yHk7iG8G3grcHXmcoj0No+WZUgYFwJD3/HyNuPi8efa9kyqdhw3TRQ/OEVuzS6/3zbXrtqKojJZA4K7X2dmZ+Usg8iYpJVRWiHncN8D/Z5dNsxUdqLZzh3bOXjLodLlo317Bqs87DON3D2ERma2G9gNsLa2lrk0Isdqe8JXn0opjJ3Hm8TK7tl3bqCsnGUpLNp+pqwcqbLVQyHwNT1HqSlma/QBwd33AHtgch5C5uKIzERTy7apwmz7jC6b0cp6BuHz6f3Sa6o2lAVxULjk6usHy0mkYZ/pjD4giCyKtPLpsjKoSTiOsm120PD3PuPs8T0O3nKodpw/beF3mYBuc+2shoAUMMopIIiMQF3LdohspvGY/YGNQ5vW/NeVadcVe7nvgSNHg8LZp59S2gMIvYmuK4Rm1aJXRd9P7mWnHwCeCZxqZrcBr3X3d+Ysk8i00tTM4RCZISqppuygQZqzqO8zrrn0gmOGmtJnxb2IoTeQDR0wNOlcL/cqo1/L+XyRsanbXdxmKWdZXp+4Ii8b829TpjgIpGWM5xrCgTl9ewqSl4aMRAaWY2IzVPLTPDtt+TftWQhBKJ1jmEWyubH2OJaNAoLIiJQlbqs7hKZuUjlVliNoiCWx8W7qeC/BNM9ShZ2HAoLIjLTZCTyvZ8eBouy0sXiSOP1M21b1CccvTnWiQFNucf4FRWZkTK3RrkMafa5P5xK6JJmrE69AuvLi83oNz2jSNy8FBJE5mcUQSjpE1LRPId59HIJAaNmHoamgbuinbNVR2z0Idd9liMAk/SkgyMoaojU6qxZs1/tVLT+tE+8vCPqsQgqfq9qd3OW7xKk6cp3XvMoUEGSpjWnIYYghlLChLIh/DuW7oi+5+vpNE9NxS7wqV1HZCqP0c/Fpa22CSNMEeDqHoUNx5k8BQVbWEEs05z3W3TfzaJAO91zwhms3vR/2EZQJu53L3k9zLg3xe9BZyfOngCBLacyTk13zC4UW/r33T4Z4QoUcdhHH16Vj++nqobodx0FVRRwmjcvmGdp+F2g/HzKGf6tVo4AgK69PxTOvDU6XXH09B2/ZPJTS92yE0MI/sHGI57zx2sqkdTffeZhdV+zdlFAvzYUU9hzE5YzLNkQ6bpk/BQRZStOM18+ycu9z/xOO37ophUV8v9BbSHsCaeqINA3GNCeVxWXqo+r7K2Dkp4AgMoVZB480CMRCHqF0Q1jZSqGqfEZxYrpwPGfoGQRlgaVNkr2qlUfpZ8q+t4JDHgoIstTaVDpN8w1VKanrlkaWDZlMM58RegLhQPogPpw+ZFVNyxTKEK9EOrBxiJ07th+95zlr2zvtAWj6nchiUkAQ6WmoHb5l4lU79z1wZFMlG3oG8RBQOsZfFXTSoadz1rYf/XvaM6gqEzTvdyjLxZS+N3TQlOkpIMjK6NoTqNqhCw+1yu+9/8imw2bqntF1PiOuVNMKNiSRS4do0gnnrsnyysoRl7lu13IcnGQxKSCIdJTuBejaU5i29ZueaRzmEcoO5KlKBzHLTV/x8tS28wRKSz0OCgiyMpoqnarX6War9FjJeFK37hlt00uE58XCOH96n1D5Vm1YK9ufkJah7v2u8yrxM2XxKCDIQhlDCzJUwqFF3nYSddpx8niJaXyf2Bajcv9A2F+QnsIW9yCGqsjrvlPf92T2FBBk5bStdOo2W0F9moZphkBChR56IuesbW/d6n7Qj813VCdedhpvOIuDSln5676zLC4FBFkIXVrXTcMZ04jH3psq3XQZZtvUDenP0/X8Ybnozh3bjw4XVWUtDcEs3WNQleAuLGPVkM9qUkAQqVBWge+6Yu/RirQuuMTDNfEyUGiXYjpcHyrouknrsgnlVFlFH6+MOnjLoU33CdTqXy3mPuUe9jlaX1/3ffv25S6GzEnXln7aiwgt6HjMPLSuu2Y2TSvM9D5lSebiyrXq8Jm6cgObnle2cid9Vro57TlvnGQzDT2IsvvEPZk4cd5QxjDvs+rMbL+7rzddpx6CLKR5VjKhYo1X8qQJ5+pa8GG4BtoNeaXDQHElXyW+5sDGIQ5sHGLXFXuPuUeZuAxlK5yqqKJfPgoIMjp18wV1SzfbLonsU5Z0jD4sNa2aeIaH0kyEz1YFjTjFdXxiWNPwTdo72GL1gSnsl0iHseIJ66oJ9D60+3jxKCDIQpm2kum71DMVKs60pxCv2om1PdgmnQyu2lVcNiyVZiBt2sncN++QKvrllTUgmNmFwJ8AxwH/xd0vz1keGYe+Szbrrk/Hy9vu0o1b7HDsWQShpxB+fsLxW48Z7okzkjalljhx29baVn7V6p84cFQN+8RzKVtscxK8WewU7nJPBZVxyBYQzOw44G3ALwK3AZ83s4+7+5dzlUnGb9pgEefbKQsKaYUbV6CxdBlnnIQu/kwIDnEgCEM7Za3z9Gd1ZwfEqSvia+smr+OT10K5ux5mrzQTyytnD+Fc4CZ3/xqAmX0QeC6ggLDipq3s69bzB2EpZ9ueQhiOSXsK4Vnxe+nn6k44q1qlFFY1VX3HXVfsPWaPwViG01JtegYafhqHnAHhMcCt0evbgKelF5nZbmA3wNra2nxKJqPXtcIoWw4at6TLln7CQ0M+TcMxVecCpAnmQuu8bH9CrMuJZG32KsRlDAEnnWDuOqegcw+WT7Z9CGb2POACd39x8fqFwLnu/ntVn9E+hOVWVSn3aeHW7QauWs+fPrts+Wea5rpqP0RdQIjF37Fu/0LVJrn4nuFeQdXvrWxfRRhKSstU9fn0d1V3fRvqGcxW230IW2pu8EkzO2vIQiVuA86MXp8B3D7D54kcrXDS1nT4+RbbvHwzXjHUphUeL9+MewDXXHoB11x6ATt3bN/0p00G1DC0lTr79FM4+/RTNt2rqUJNd0CfuG3r0Qlmkboho3cDnzKz9wBXuPsPBn7254HHm9ljgW8Azwd+feBnyAKZdrIytLLjvD9l92pb+ZUdbB/uOYuWbNnO4XTpaJt02l3Llh620/T5Wa5IkrwqA4K7f9jMrgVeA+wzs/cCD0bvXzXNg939iJn9LrCXybLTd7n7l6a5p0iqLLVz2aHv6VBN6A2kZxi3nYyuGkZJh6nKgldcpqbgVva6jlYISZ2mSeUfAPcC24CTiQLCENz9k8Anh7ynLL5pKqmyVnPVkEsb8UqfuvMPhswOmvYUYrM86axPj0yWS2VAKDaNXQV8HHiKu39/bqUSKdGlVZtO4taN/ZdVwGWTzmEitir/UPy8dIK36bnx6/j9uuGjaagylzJ1PYTLgOdpGEcWSdW+gHvvP3LM5rKu90xTSIeVPmXnEUwrHd46eMuhTXsPhszXJBLUzSH83DwLIlKlz+alsmWcVQfRx58pe24XZfet07YSj4euFoGC02JarP+XiVSoChpxTyFeNRQyf7ZRtZs4brHDQzuGh1jCWTWcVNcz0G5fmZYCgozeECtjuozB161Gavucvs9Y9Ep8Wb/XqlBAkLlqW0FU7SiuUtWKT++RTv6G3b5tTgir2gPQpZxdxXsrqjaeaSmpDEUBQRZGmzQR6c/DjuFQ8YfhnLINZ+lzZlnBdn3GLJebDknBabEpIKywef5H23YoIZ0IDpV5l55CWU6gcCZBuEeaB2ia3dFl33OI32nblN1VZRHpSgFBFlJd4CiTLhlto0sF23TfvruMp0nZndPYyyflFBBW0Cwn/poqvjbPik8p6zM2H/YHpGkn4nuGOYNpvnv82TQz6hC/07LejpLQySwpIMjcxOP9Q1SY8TLSNE10utwUygNM2+ylVeUMLfaqQ3CmDb5Vk+Mis6CAsIKGmvgrO3+gqeKrOmimLMd+m7MHqsoUB4qDt0z+/qCzKSX1lRefN8gh8+H+QdiPMFSqiSFTVojUUUCQTrqsigmVbVqBlqVg6Kpp41YqPcqyKZNol5b9Ccdv3TRcFO8qbht8m95Xz0DmQQFhhfWtZOoq0zbDG+FQ+HAvoHRMv83ZBlXfKZQjrC5KU1ZMMxYfyrbFJpX/NZdeoGEdWQoKCNJK2RLIsmvSc3qhfBNXvA9g2sq5TDzJWza+Hyr1dO4hvW8Y2mqq5JuGdZp6BtrZK2OggCC9hEnb9LzfqlPG4s/FWULTyjjdmRvuGwJQmzmFODtpELfmp5Gu+rn3/iOd9kmIjJkCgrRSlRqi6Zq698smeVNxMCg76KZPizptlYeflSWMm3auo4l29sqYKCBIb0MtGQ3Khk/CoTShYj54y6FNa/3jyev4HqEHEz53zlrzgTVtzCuHkUgO5j7wyR4ztL6+7vv27ctdDOmgqfUfv5cGhLKNZbETt20tPaEsDhJ1lXabVnlTig0FBFkEZrbf3debrtsyj8LI4gjj9/MSDwOFbJ6hck/nAIITt21li22ejL75zsPcfOdhDmwcOpriAWa3s7fvHoYm8/79i8Q0ZCQzUTb8E69MqjqfIL0uXiaaDg+lny/L+xM+W6ZLWm2RVaCAIMD8lz9W7WVId/2GIaGm1UVVeX/anoqWajPUVbUPo+pzbZ6n5aeSkwKCzEQ6+Zqmo4ZjK/8y8VBR2SqgsiWraU+haSWTiEwoIAjQfSPWEELKh/icgrSn0HV1UFlQ6KLNOQlVS3C77qpOy93ncyJDUkBYUVUrfNocGNOl0kpb9XE6iWsuveBoBVy2NyB+nc411FW8VT2JdMObiGymgCBzk6aTiJPcHdg4dMxKIihPOzELaUCsGqqKNfUa+lCwkpyyBAQzex7wOuCngXPdXZsL5qRp8rLLZGrXlA115xmH1BLhurryxInl2uYYKsuwWvfZuv0PIssq1z6Eg8Au4LpMzx+9ZV6PvnPH9qN7DUJLPPQa4mMihxJ6JnVCcAg9kyBMSLctT9hL0dcy/7vL+GXpIbj7jQBm1nSpDKxtvqGqz4YMoaHyDvcacqgj3WsQ9Ekslya6CxV+m9TcQbr/QWRZjX4Owcx2A7sB1tbWMpdm9sa+Hr1NazuVfqfQOwiZR8PEcqh0h/quaQCB8go/KNvLUJUee2hj/3eX1TCzgGBmnwFOL3nrMnf/WNv7uPseYA9MchkNVLyV17eiSVvKbe7TZgnofQ8caWyJ90ksV3fuctUz4mWrylUkq2RmAcHdz5/VvZdZjv0AbaQt2HSsvU5ZEImXgIYTx8J7Q2kKIHXDZvNuoWsfgozB6IeMVk3ZfoBFzarZJoj0GSrp83vo+vtbtN+1yBByLTv9ZeA/AKcB15rZF919uqOsZGbqDsWpuzYVcgvFO3q79DT6KktdDeMcrx9DGWR15Vpl9FHgozmePXZlx0Z2zcczywouLlPTs9qcslY3JzGmilpkFWjISCpVtaTbnnYGzcdopq/nQeP1IuUUEEaqT0t5HkMh6TkGbSr0NnmEyr7vGId0RJaZAsKK6jJ5Gw9fNVXobSvvMVTuYyiDyJgoICyAPimU47X/s8jw2fZ+XVr7XfIqicjwFBBWTNfhmHgZbJuJbVXeIotLAWHJpBV4+vN5m0cWVREZhgLCiuk71r+Kwzer+J1ltSkgLJmqVT+5K7UxZFEVkXoKCCMyzxZp12c0zTEs04YyLXuVVaWAsKS6Vl65K71VHpoSGQsFhIFMU5HVtUjHWkGWlTkk4VvklvUi/O5FZkUBYcl0rcTGNjyiylckHwWEKQ1RoVYlgYuzguauqFN1QzxjSMLX1dgCo0gOCggZ9Klsmk4d61uhrcLY/TJ/N5EhKSBMacgKNf5sjiygfdSlva4yxtb4KgRGkSYKCHPUpyJse/7AtBXaMlaAYww8ImOmgDCQWVUy4QCZ9PyBRTbm1viYyiIyb+buucvQ2vr6uu/bty93MabWpyLskiV0Ucyr3Iv6+xEZipntd/f1puvUQ5BsVEGLjIt6CBUWoVWZjpHv3LEdGHeZRWT+2vYQtsyjMFIu7DUQERkDDRklFmFlik4WE5FZUEDIYBGCjoisHgWExJhb3YseSBatvCKrRgEhgzEHHRFZXVpllFHfU80WLZBoNZRIXqNeZWRmV5rZV8zsgJl91MwemaMcuV158XmqFEVkNLL0EMzs2cBfufsRM3szgLu/oulzy9pDmEXLeYy9iDGWSWQVjLqH4O6fcvcjxcsbgDNylENERB6SfQ7BzP478CF3f1/F+7uB3QBra2tP3djYmGfx5mLIlrPG60UklT2XkZl9Bji95K3L3P1jxTWXAUeA91fdx933AHtgMmQ0VPk0fCEistnMAoK7n1/3vpm9CLgIeJbn7qZkNmRQqlrSqgAoIk2y7EMwswuBVwDPcPfvz/PZi765S0RkVnJtTHsrsA34tJkB3ODuL81UlsGMKbikPQMFQBFpkiUguPvjcjwXNIQiIlJFqSsGMOZWuAKgiLS1sgFBFaOIyGbZ9yF0MfadymqFi8gYjXqnsoiIjM/KDhnNgnoGIrLI1EMQERFAAWHhXXL19cecqyAi0ocCgoiIAJpDyGbaFUlj3vsgIotJPQQREQG0D2Hu2pxX0KW1r56BiDTRPgQREelEPYRMylr2Ou1MRGZBPYSElmeKiNTTKqNMylr9ykwqIjktfUDQ8kwRkXaWPiAsIgUrEclh6QOChmFERNpZmUllERGpt/Q9hEA9AxGReuohiIgIoIAgIiIFBQQREQEUEEZDO6lFJDcFBBERATKtMjKzNwDPBR4E7gJ+w91vz1GW3LSTWkTGIlcP4Up33+nuTwI+AbwmUzlERKSQpYfg7oejlycCi5ODe2DaSS0iY5FtY5qZvQm4GPgH4J/mKoeIiEzM7IAcM/sMcHrJW5e5+8ei614JPNzdX1txn93AboC1tbWnbmxszKK4IiJLq+0BOdlPTDOzHcC17n5O07XLdGKaiMi8jPrENDN7fPTyl4Cv5CiHiIg8JNccwuVm9gQmy043gJdmKoeIiBRyrTL6lRzPFRGRatqpLCIigAIqIniCAAADRklEQVSCiIgUFBBERAQYwbLTLszsW0wmoZfFqcC3cxdijvR9l5u+73jtcPfTmi5aqICwbMxsX5u1wctC33e56fsuPg0ZiYgIoIAgIiIFBYS89uQuwJzp+y43fd8FpzkEEREB1EMQEZGCAoKIiAAKCNmZ2ZVm9hUzO2BmHzWzR+Yu0yyZ2fPM7Etm9qCZLdWSvcDMLjSzr5rZTWb273KXZ9bM7F1mdpeZHcxdllkzszPN7H+a2Y3F/49/P3eZhqSAkN+ngXPcfSfwf4BXZi7PrB0EdgHX5S7ILJjZccDbgOcATwR+zcyemLdUM/du4MLchZiTI8DL3f2ngacDv7NM/74KCJm5+6fc/Ujx8gbgjJzlmTV3v9Hdv5q7HDN0LnCTu3/N3R8APgg8N3OZZsrdrwMO5S7HPLj7He7+heLv3wNuBB6Tt1TDUUAYl98E/iJ3IWQqjwFujV7fxhJVGPIQMzsLeDLwubwlGU6uA3JWSpvzpc3sMibd0ffPs2yz0PY87SVlJT/T2u4lY2YnAR8BXubuh3OXZygKCHPg7ufXvW9mLwIuAp7lS7AxpOn7LrnbgDOj12cAt2cqi8yAmT2MSTB4v7tfk7s8Q9KQUWZmdiHwCuCX3P37ucsjU/s88Hgze6yZHQ88H/h45jLJQMzMgHcCN7r7VbnLMzQFhPzeCpwMfNrMvmhm/yl3gWbJzH7ZzG4DzgOuNbO9ucs0pGKBwO8Ce5lMOH7Y3b+Ut1SzZWYfAK4HnmBmt5nZv8ldphn6GeCFwC8U/71+0cz+We5CDUWpK0REBFAPQURECgoIIiICKCCIiEhBAUFERAAFBBERKSggiPRUZL78ezPbXrx+VPF6R+6yifShgCDSk7vfCrwduLz40eXAHnffyFcqkf60D0FkCkUag/3Au4CXAE8uspyKLBzlMhKZgrv/wMwuAf4SeLaCgSwyDRmJTO85wB3AObkLIjINBQSRKZjZk4BfZHJ61h+Y2Y9nLpJIbwoIIj0VmS/fziQn/i3AlcC/z1sqkf4UEET6ewlwi7t/unj9H4GfMrNnZCyTSG9aZSQiIoB6CCIiUlBAEBERQAFBREQKCggiIgIoIIiISEEBQUREAAUEEREp/H+CJu4anYEAdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[0,:,0],Y[0,:,0],label=\"data\",marker=\"+\",color=\"steelblue\")\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title('Plot of data samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the MINEE model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from model.minee import MINEE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 100        # batch size of data sample\n",
    "ref_batch_factor = 300  # batch size expansion factor for reference sample\n",
    "lr = 1e-5               # learning rate\n",
    "\n",
    "minee_list = []\n",
    "for i in range(rep):\n",
    "    minee_list.append(MINEE(torch.Tensor(X[i]),torch.Tensor(Y[i]),\n",
    "                            batch_size=batch_size,ref_batch_factor=ref_batch_factor,lr=lr))\n",
    "dXY_list = np.zeros((rep,0))\n",
    "dX_list = np.zeros((rep,0))\n",
    "dY_list = np.zeros((rep,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Previous results loaded.\n"
     ]
    }
   ],
   "source": [
    "load_available = True # set to False to prevent loading previous results\n",
    "if load_available and os.path.exists(chkpt_name):\n",
    "    checkpoint = torch.load(\n",
    "        chkpt_name, map_location='cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    dXY_list = checkpoint['dXY_list']\n",
    "    dX_list = checkpoint['dX_list']\n",
    "    dY_list = checkpoint['dY_list']\n",
    "    minee_state_list = checkpoint['minee_state_list']\n",
    "    for i in range(rep):\n",
    "        minee_list[i].load_state_dict(minee_state_list[i])\n",
    "    print('Previous results loaded.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuously train the model. The following can be executed repeatedly and after loading previous results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VFXawPHfk05ICC30jjQRUYiCUkRFRVHRd9UXQde1r67r666r4trL7uIW2+q6i713XXsBRbEAElSagLTQS0IJgZA65/3j3oTJ9Japz/fzyScz955777kzyTNnThVjDEoppRJfWqwzoJRSKjI0oCulVJLQgK6UUklCA7pSSiUJDehKKZUkNKArpVSS0ICeBETkCxG5LMZ5GCAiP4hIhYhcG0D6O0XkBftxDxHZJyLpzZ/T5CUiU0Xk01jnQ8WOBvQEISIlInLADnzbReRpEckL8hy9RMSISEYzZPFG4AtjTL4x5uFgDjTGbDDG5Blj6pshX0nJ03tpjHnRGHNyM10v5oUG5Z8G9MRyhjEmDxgGHAXcGuP8OOsJLIt1Jpw10weXUnFLA3oCMsZsBj4CDnPdJyJpInKriKwXkR0i8pyIFNi759i/99gl/WNE5BAR+VJEykWkTERe9XZdETlTRJaJyB67xDbI3v45cDzwiH3e/h6O7W1fp0JEZgLtnfY1ljZFZLKIFLsc+zsRedd+nC0ifxeRDfY3lX+LSAt73zgR2SQiN4nINuBpe/uNIrJVRLaIyGX2tQ4J4nzX26/lVhG52ClfLUTkH/ZrXS4iXzsdO1JEvrVfq0UiMs7H69pFRN4UkVIRWedcZSUiR4tIsYjstfN3v4/38lci8rXTsUZErhaRVfbrfo+I9BWRufb5XhORLDttGxF5387DbvtxN3vfn4AxTu/vI/b2gSIyU0R2ichKETnP6dqnichP9nU3i8gfvN2/iiBjjP4kwA9QAoy3H3fHKg3fYz//ArjMfnwJsBroA+QBbwHP2/t6AQbIcDrvy8AtWB/uOcBoL9fvD+wHTgIysapYVgNZrnnwcvxc4H4gGxgLVAAvuOYLyLX39XM6dgEw2X78IPAu0BbIB94D/mLvGwfUAffZ12kBTAC2AYPtcz9vX+uQIM53t33PpwGVQBt7/6P2fXcF0oFj7et2BXba6dPs12wnUOjhdUkDFgK3A1n2+7YWOMXpdbvQfpwHjPTxXv4K+NrpubHvrZV9/9XAZ/Y1CoCfgIvstO2AX9ivUT7wOvBfp3M1eX+BlsBG4GL7fRsGlAGD7f1bgTH24zbAsFj/D6XCT8wzoD8BvlFWQN8H7AHWA/8CWtj7Gv/Z7H/Yq52OGwDU2v90noLAc8AMoJuf698GvOb0PA3YDIxzzYOHY3vYgbGl07aX8BDQ7ecvALfbj/thBfhcQLA+VPo6necYYJ39eBxQA+Q47X8KO0Dbzw+xr3VIgOc74PJ67QBG2vd/ABjq4X5vwv4Qddr2SUPwdNk+Atjgsu1m4Gn78RzgLqC9SxpP7+WvcA/oo5yeLwRucnr+D+BBL+/ZEcBup+dN3l/gf4GvXI75D3CH/XgDcCXQKtb/O6n0o1UuieUsY0xrY0xPY8zVxpgDHtJ0wQr4DdZjBfOOXs55I1Zg+86uTrnES7om5zXGOLBKaF0DyHcXrOCw3yVf3rwEnG8/noJVUqwECrEC+0K7KmMP8LG9vUGpMabK5dobnZ47Pw7kfDuNMXVOzyuxSsrtsb7RrPGQ/57AuQ3ntM87GujsJW0Xl7R/5OD7dSnWt6MVIrJARE73cA5ftjs9PuDheR6AiOSKyH/s6qO9WB8krcV7z6OewAiXfE8FOtn7f4H1DWW9XdV2TJD5ViHQRqPkswXrn61BQ+l4Ox6CrzFmG3A5gIiMBmaJyBxjzGoP5x3S8EREBKvqZ3MAedoKtBGRlk5BvQdWCdKTT4H2InIEVmD/nb29DCsIDTZWO4InrufcCnRzet7d6XEg5/OmDKgC+gKLXPZtxCqhXx7AeTZifSPo52mnMWYVcL6IpAH/A7whIu3w/tqF6nqsb3MjjDHb7Nf+B6wPezxcbyPwpTHmJC/5XgBMEpFM4BrgNZq+9qoZaAk9+bwM/E6sRsg84M/Aq3YpsxRwYNWhAiAi5zY0fgG7sf5xPXUffA2YKCIn2v+k12PVyX7rL0PGmPVAMXCXiGTZHxxn+EhfB7wB/A2rbnumvd0BPA48ICId7Px3FZFTfFz+NeBiERkkIrlYddUN1wnlfM7HPgXcbzdqptsNk9lYVUZniMgp9vYcu4G1m4dTfQfsFasht4Wd/jAROcrOzwUiUmhfb499TD0e3ssw5WN9uO0RkbbAHS77t7tc632gv4hcKCKZ9s9R9uucJVaf+AJjTC2wF89/UyrCNKAnn6ewGv7mAOuwSpG/BbCrLf4EfGN/TR6J1f1xvojsw2pA+z9jzDrXkxpjVgIXAP/EKp2egdWNsibAfE3Bqi/ehRUsnvOT/iVgPPC6S5XHTViNsfPsqoFZWCVLj4wxHwEPA7Pt4+bau6pDOZ+LPwBLsBptd2E1xqYZYzYCk7CqTkqxSrM34OH/zVh978/AqrNeh/XaPoHVaAlWo+4y+/15CKtxuMrLexmOB7EakcuAeVhVT84eAs6xe8A8bIypAE4GJmN9e9vGwcZogAuBEvs1/TXW345qZmKMLnChUodYXS2XAtkuHxRKJTwtoaukJyJn29UAbbBKke9pMFfJSAO6SgVXYlV9rMGqy70qttlRqnlolYtSSiUJLaErpVSS8NsPXUSeAk4HdhhjDnPZ9wesrmWFxpgyf+dq37696dWrV4hZVUqp1LRw4cIyY0yhv3SBDCx6BngEl25mItIda46KDYFmqlevXhQXF/tPqJRSqpGI+BpZ3chvlYsxZg5WH1tXD2ANG9dKeKWUigMh1aGLyJnAZmOM65BnT2mvEGv6z+LS0tJQLqeUUioAQQd0e/j0LTgNofbFGDPDGFNkjCkqLPRbBaSUUipEoZTQ+wK9gUUiUoI18dH3ItLJ51FKKaWaVdCzLRpjlgAdGp7bQb0okF4uSimlmo/fErqIvIw1odEAsZbjurT5s6WUUipYfkvoxpjz/ezvFbHcKKWUCpmOFFVKqWZgjOG14o3U1Dmidk0N6Eop1QzeW7yVG99YzCOzXRf/aj4a0JVSqhmUH6gFYOe+aj8pI0cDulJKJQkN6EoplSQ0oCulVDOK5mRXGtCVUipEm3ZX4m2RIIlyXkADulJKAbCnsoaft1c0Pi8u2cU7P272mn7QbR8z+r7ZPPNtSRRyF5igh/4rpVQymvToN6zfWUnJ9IkAnPPvudb2I7qyt6qW79buYvyhHek17QN6tsvlQG09AAtKdnF4t9ZsK69i4uGdY5Z/0BK6UirOzfxpO6ucSs7BWrFtL6t3+D9+/c5Kj9t37a/huld+5LLnitm854BbWmPgF499y29e+j7kPEaKltCVUnHt8uesVc4aSs7BmvDgV2Edf8ET86myS+MNv4PhpYq9WWgJXakUcdmzC/jD637XpEkq4/42221beWUtj85e7bUx09VPW/eytmy/1/3Op9lXXccTX63FGIPEoFVUS+hKpYhZy3cA8Pdzh8Y4J5Hz25d/oN7h4F9Thzdu27S7kvycTApaZFLioRrl1neW8t6iLRzerYAx/dwX3SnbV037vOyQ8jPs7pnU1Dtol5cV1TlcGmgJXSkVtypr6po8f+W7DRzzl88an7+3aAsfLtnWJM3o+2Zz2kNfeT3nviprSH5t/cGAe+t/lzQ+Lrp3VsCldwDj1NO8xj7n715dxF3v/RTwOSJFA7pSKm795sWmDY3T3lrC1vIqt3Srtlfw7y/XND5vaLx09sRXaz1eo7qunhfmbWiy7ds1Oz2mDa0+PHqV6BrQlVJxa/bKwBaWP+mBOUz/aAV19d6rOe79YHmT58u3Wj1f7n1/uVva6rrAGz8ra3yn3VcdfENqqDSgK5WCNu7y3EXvhw27OX/GvJDqfx0Ow/uLt+BweC6RVtXWc8BP8Ht30RbmeikdB6LSqReKr+AO8LdPVgKwba97id+b8fd/6bbtq1W+V998b9EWek37IOBrhEMDulJJbvnWvbz9w6bG558t386Yv87m46Xb3NLe9OZi5q7dyTofvTqK7p3Fb1/+wW37i/PXc81LP/DKgo1u+5ZsKmfoXZ8y6PaPvZ73m9VlXPvyD5z/+Dw27fb8gdOgbF81b32/yW374Xd+2vh4X3Wd2/5ARbOrYSRpLxelktypLg2Ey7bstX+XM+GwTo3bN+ys5Oft+7yep7KmjkNv/wSwSp3/PP/IJvt3VFjzfj83t4Tzj+7Oif/4ktH92vPc3PUB5XPqE/MP5mVXJd3a5HpN+8p3G9waQ12Jh36Dnr49eAreGtCVUgltrIc+284+8hNAG6zYVkG/Wz6izmF89t/258X5TT8InBs6/dSmAJDmoR94nz9+6CGle/R+ZcEGD+nin98qFxF5SkR2iMhSp21/E5EVIrJYRN4WkdbNm02lVDSd8uAcXvmuaVC7PohBSXVe6tHBqksvP1DLs9+W8NGSrQDUu6Sf8vh8bnl7aZNto6Z/3vh4+da9/LR1r888pAUwsqfXtA/YU1nrtr2hz36iCaSE/gzwCPCc07aZwM3GmDoRuQ+4Gbgp8tlTKjGVV9byybJtnHdU91hnJWTT3lrC5KN7BJze21work596CufdfSB+HiZ/28LgY7ULF6/O6y8xBO/JXRjzBxgl8u2T40xDS0O84BuzZA3pRLW71/7kRvfXMxyP6VIV09+vS7kYPfz9oqgBsSUH6hl6ebygNJ+tny727Z/fdF0+Py7i7YEdK5wg3mgfLUHxEIwXSFDFYleLpcAH3nbKSJXiEixiBSXlgbWp1SpWPt5e0XAwc6TUnth4Fe+28CijXvc91dUu3UNrKqt5573f+Kcx771e/5V2yvoNe0DvlltdZn7dk0ZJz8whxfm+6/7vX/mzwA8N3c9p//za7/pAb5Z7d6V8K8fr6T3zR9y6TMLAjpHtJ316DexzkITW/cE3j0yVGEFdBG5BagDXvSWxhgzwxhTZIwpKix0nzdBqXh08gNzAgp2g277mJveWOx1/7Nz1zPJQ2A56k+zuNap61+9wzDwNqtLn3N3u9krdjTWMzubv8760vyBva+kzKru+HFD0w+PtaX+S6m/evo7n/3DF23cw75q93rmBp+t2MHo+z73ul9FT8gBXUQuAk4HpppgvucplUQO1NbzarF7v2tX93+6kpNcBqU01ANv3nOgybwiBqipc1BdV8/FzyzgKqfh77X1DnpN+4Db3mnaYNjQb/tNu2/2xl2VlFZUc8I/3AfCuPpiZSmvebmHh2atYtKj3/BasXuf76bXdx9qr5oq2dn8VU0hdVsUkQlYjaDHGWMCawlRKsE98806WudmsX5nJYd3K+D4gR3c0hhj2LjrgNt6kg9/vtrjOT9dto0rnl/Io1OGNW6rqXPQ/9aPyEp3L2+V2n29XYtQu516atzw+iJeX7iJq8f1DfDO4I53l3nc/sCsnwM+h/Lt9eJNjBvg/jcTSX4Duoi8DIwD2ovIJuAOrF4t2cBMu/P+PGPMr5sxn0rFxNbyA/zny7Xcdvqh3Okye57zggkOhzX/9RsLN3GDjyoYgHveP3ieJXY9/RIP9fU1gXS2xqp7f9mpi+HrC63S9MIk6r2RFKIwP7rfgG6MOd/D5iebIS9KxZ1j/mLVDZ8yuJPbPuf5Ofr88UPOPrIruVnpPs9XXlnLk1+vc9vuPFOgL7v21zR5/tL8DdR6mXdF60HjSzTWu9C5XJQCXv5uA72mfRDW/B9v/+B9hfgGQ+/+tMnzx74ILJA3lLbfWOhel/26h20A363b5XG7io1ABjqFfY1mv4JSCeDxOdZc2dv3VrmNWgRYWxZYn+YXA+g26MzXiEpnv3jsWz5Zto33FwfW11vFH09TEUSazuWiUs66sv1s31vFyD7t3Pbt3FfDQ7NWuW2fMcfz4gjRdOXzC2OdBRUGT5OFRZoGdBVV28qrmL9uJ5OO6BqzPBz/9y8Aq1Fz1k/bOaJH68ZJpM77z1yPxwQ6rF0pb1Zsq2j2a2hAV1E15fF5rC3bz8mHdqKFjwbEiqpa1pbuZ2h33/O+rSvbT+eCHIpLdrNtbxXnDD84C8X//OsbCvOzqas37K2q5fVfH9vk2PIDtVz2XDGDOrcK76aUCoC3RUUiSQO6iqqG9SCNnz4Ylz9XzLy1u1h57wSyMzwH/sqaOo7/+xdMPLwzHyy2Rkw6B/TvN7gPuXc29C6rgTLY+VaUCkU4De6B0kZRFRP+xhb/YAfjb1aXMXuF56lMq2sdjWka3PnuMq89R47/+xfMXpmY06IqFQgN6Cpq9lbVcqDWfc6Q8gO1nPPYt42rshtjqLb7Vl/yTDEX25M/7auuo9e0D3hhnrXwgafPhGe+LeG+j1d4vP66sv1c/HR8TiSlVCRoQFdRs77Mcx3iw5+tonj97sZV2b315NtuL+Z763+X0mvaB/z5Qyu9p74DG7QRU6UgDegqapzrzZ1jtnO/b2+B2HmR4waeBtk08LecmlLJSBtFVcw5T9Y59m+z+enuU9zS/O5V78uf7fawhJhSqUhL6CpqnKeILXFatca1hqVhZXmlVHA0oKuQOBwmoCW1vl1TRl29g9+/9iO/eOzgoJ3T//k1364pwxiDQ6fTVyoiNKCrkPzx7SUMuNVaYafKQ88VsCaHmvL4fB76bBVvfe8+cdWUx+fz6xcW+u3CqJQKjAZ0FZJXFlgr3Dz19ToG3vaxx8E5K7ZZ2/7pZXEHgE+WbddpXpWKEA3oqlFNncPnqvFVtfVMeXweK53mpLjbXqxh2RYreM+Ys4avV5Xx/Lz13P6O51VwXOkKhkpFhkTzn6moqMgUFxdH7XoqcFvLD3DMXz7n3rMO44KRPT2m+WZ1GVOfmE9WRprbivVKKd9EYN1fJvpP6PFYWWiMKfKXTkvoKezhz1Yx9q+z2bS7snHV+HcXNZ1v2+EwvFa8kV7TPmDqE/MBNJgrFYJorFik/dBT2P0zrQWAR983m5cvHwlYAXzzngM8Pmctw3q2YfqHy9liT6illApdNFYs0oCeIkorqinZuZ/qWgdDuhbwt0+bzneycbdVQncYw6jp1jqaz3xbEu1sKpW0ohDPNaCninF/m83+Gqt74aQjuvDOj02rVm60V6r3N+WsUio0EoVKF7916CLylIjsEJGlTtvaishMEVll/27TvNlU4WoI5nBw2lmlVPREo4QeSKPoM8AEl23TgM+MMf2Az+znKkGUH9C5T5SKtrgI6MaYOcAul82TgGftx88CZ0U4XypC6h2G9Tv3N9k2d+3OGOVGqdQVjSqXUOvQOxpjtgIYY7aKSAdvCUXkCuAKgB49eoR4ORWs0opqrn99EUs3l7Nrf02ss6NUyouLEnq4jDEzjDFFxpiiwsLC5r5cSpr503ZWbT84etMYw29e+p45P5dqMFcqTsRzP/TtItLZLp13BnShxiiprqsnIy2N9LSDfx6XP2eNvj2qVxtOPaxz43B8pVT8SEuL3yqXd4GLgOn273ciliPlU8MMhyXT3YcQLyjZzYKS3dHOklIqANEYWBRIt8WXgbnAABHZJCKXYgXyk0RkFXCS/VwppZQXBS0ym/0afkvoxpjzvew6McJ5UX44nNbe/GTZNq58fiEjereNYY6UUoGafHT3Zr+GTs4Vx3bvr2mybNvt7zaO7eLK5xcCMH+da49SpVQ80rlcUpjDYTjynpmcfWRXKqpqmbOqTGc5VEr5pAE9Djkcht2VVnfDt39wX7pNKaU80YAeJ14v3shJh3aksqaeY+3ZDpVSySOe+6GrCFq9o4Ib3ljMcf0LmTJCR9MqpUKjAT2GjDGU7avh5e+sBZfL9lVT79D1NZVSodGAHkNnPPI1SzfvbXxeUVXHoo06H7lSySgp5nJR3jkHc4ANuyr5z5y1McqNUirRaQk9BhwOw32frPCfUCmlgqABPcoqa+o457G5/LR1r//ESqmkEc/zoasQnfzAHDbtPhDrbCilkpDWoUeRw2E0mCuVogzN34NNA3qUGGM47eGvYp0NpVQS0yqXKPhu3S7O+8/cWGdDKZXktIQeBRrMlVLRoAG9md3xzlL/iZRSSU97uSSw6rp63v1xC8/OXR/rrCilUoQG9GZQWVPHobd/EutsKKVSjFa5NIOT7p8T6ywopVKQBvRmsHmP9jVXSkVfWAFdRH4nIstEZKmIvCwiOZHKWKLaX10X6ywopVJUyAFdRLoC1wJFxpjDgHRgcqQylohWbNvL4DsiU3f+0uUjInKeYPz2hEOifk2lUkUiTJ+bAbQQkQwgF9gSfpYS14QHIzcS9Ni+7bn+pP5Nto3s05Zvp53g9Zg+7Vs2Ps7PDr69+/qTBwR9jFIqMCYKa9eEHNCNMZuBvwMbgK1AuTHmU9d0InKFiBSLSHFpaWnoOY1zP2+vCCjdFWP7NHn+5Q3j+OS6sQBceVzTfa6f6E/96ii6tG7h9dz5LTIbH2ekey8OHD+gkGvjsDTeva33e1NK+RdOlUsbYBLQG+gCtBSRC1zTGWNmGGOKjDFFhYWFoec0zp38QGA9Wy4c2bPx8dc3HU/Pdi0Z0CmfkukTufnUQdx2+qE8MuVIACYc1rnJsblZnkvdR/Vqw60TBzHjwuGN2xoKAw3ncvbERUeRmR7cW3/XmYODSu/sj6cNDChdejS+kyqVxMKpchkPrDPGlBpjaoG3gGMjk63E8s3qMr9p1v75NFb96VQ6tMqmfV4W/75gON3a5Lqlu3R0b04/vAsAh3TIo2T6RLc04wd1aFKlMqhzKy4b04eOrXL4w8n9ee+a0bTNzQJwK9H375hHepqQnWm99VNdFqWeMLgTAEO7t26y/aJje/m9R29at8gK+dhoa5GZHussKBWycAYWbQBGikgucAA4ESiOSK4SzNQn5nvcPueG4xn7t9m0zs0kLU1Is4f+Ft96UljXe+KiowDYuKuSJZvLOWFgh8Z915zQD4DnLxvBFyt3UOBUDfPMxUcxboCVtqE+LzeraQD7t13KX7KpnDMe+RqAc4d3A6B3+5asK9sfcD7POqILP2/f5zPNrN8fx/j7vwz4nM1p4a3jyc5M57pXfmDW8h2xzo5KMl3bNH+VYsgB3RgzX0TeAL4H6oAfgBmRyliimPnTdo/bp47oQY92uR5L2MG6deIgxEN1RPe2uXRv617KB+jaugVTR/Sktt7RuM1Tm4yI8NjUYbR0aUTt0c46b+/2LfnbuUMBeOeaUZzywBy2llfx0uUjmPK45w+yBg9Otqp7Xluwscn2UwZ3ZPWOfawp3c8hHfJ8niOa2uVlA/DIlGEMvO3jGOfG+iamHyzJIxoVimEN/TfG3AHcEaG8JKTLn3P/UrLkzpNp6aW+OxSXjenjP5EXmelpnDK4I58s2+61jvrUIZ3dthW0yHT7MGqVk0mPtrlsLa8KOT8A/7mwyOf+od0KOG1IZ576Zh3b91YD8NSvirjkGd9fAPsUtmRtaeDfILzJ0WoXlaB0LpcwPPzZKo/b83MyPW6PlT+fPYR+HfIZfUj7xm3h9qAKZua4YFZqWXHPBDLShIz0NK48ri+9pn0AQEEA9fDPXnw0Y/46O+BrOXto8hEhHadUPNGh/2G4f+bPbtvC6Q3SXNrlZfOHUwaQluYehIP9Gvjr4/oCMKhzPmP6tW+y7+JRvRofn3pYp2CzSfe2ueRkppMRZA+cBt1CrKNccMt4Jh3Rtcm2o3q18XlMVx/dRyNHe/2o4GhAD1G9w3Op07lbYjI6fmAHSqZPpHVuVpMGV4A7zhjcWNL989lDPB7f3q6n9uSR84eFlTdP7Qy+jO1vdaMtzHfP05Vj+/o89uXLR3LrxEFBXS9Yrq+vUv5oQA/RPz9vWt0ypGsBJdMneiwFx6NIjFr701lDmHZq0z7mk47oSsn0ibRp6V5Fcu7wbhTfOt7r+Qpy3QPYB9eO5ssbxoU9bPrGCe6jYKdNGBhyo3WPdrlhtW0EIpIDre4567CInUvFLw3oIXpwVtOA/u41o2KUkzCFESgLcjMbq2ACulQI1xrcpYCe7Vr6T+jH1eMO4fbTDw37PIkq2b85JoIojPzXRtFQOFyqW1bcMyHor/vJ5Mfbfferb+haObBTq2bLw7+mBl9d4+st6xtH3SmVCpSW0EPQ548fNnmeiN3cTh7cEYBJQ7v6Self69wsWud674VybN/2vHfN6CaNppHm2kAbrt7tW7LsrlMCTt+rnefxAOGIxhqUKrloQA+Scal8jsTAoVjoW2hNK3Bol+YrNTsb0q3A57eYiYe794V35iu0Lb97QrN0FXUdbOVLMF+n3//taD64dnTj8zd+fYzXtAt9tDkE6+2rU3JmjpSiAT1Im3brakSRtu4vp/HI+e6TiAWqRVZg35Bc56eJJOfP+e9v810FdVjXAgZ3KWh8XtSrrcd0IgdHr0bCkT18d8UM1xSXeYFUU9H4vqUBPUjOA1dce3io0IiI3zaIhv1DuxUwvGdogSnU43xp6LqYlXHwX6mthx4+/sz6/XERy1OseOuqqizRaBTVgB6Gy5u525pyZ4hc/+xg2rE/uW4sK+6Z4La9oS+7YPVNf+bio0LKi6c5bYIt0Q3qHJ3qMxW/NKAHwbl3y+hD2pOeIH3Ok0E4r/T4QR3Dvv6ATvk+G78NcEzfdo2zWUbbqYd14qP/G+M33SWjekchNypWNKAHwbl3y4xfDveRUkVLIKM1G+Z+bw7x9pHub8qC289I3b74qUADeoi8rR6kmpcxB4PoE78savbRmv40LCDym+MDH2AVqFCGNrz+62Obpa1AJQYN6AFy7q541hFdYpiT1OQc3IbZAatz65yAjm2V4/nDN5CujicO9F2F0jI7g5LpEzn7yG4e979w6Qheu9Lqltjaw9QG4QqlATYaOjjNj9O5IMdnz5/rxveLRpZSghYzAzRn1cFl5v5uL/igos9guOq4vpwyuCOHdMj3m/7WiYOYfLR7d7r3rhkd0IyJT/4qtEbOBiP6tCUzPY3ld08gJzPNZ2+ev55zOPe89xMV1XVA8JONNQgh1ldJAAAZFklEQVTlqKz0NGrqHVw5tg9l+2p48/tNQR3/W5dFxwd2bsWOitLG/Pj64LlufH+3qTSSUSTmT/JHS+gBuuip7xofhzq9qwqd86jJtDTxGsyP6N6ah536tF82pg95HgYIDelW4LYtEOcVeS6J+9MiK91vgD6vqDu/PPbgnCveVqNqDkf3tvrCj+7XPuyJ0OLFxCGdOW1I8NM4JzKNTEF6ZEroA2BU6NrnWyW8UX19D/H/729GceZQ71Vir14xkpcuHxFyPv56ztCojA4e0689Z9ijZ3sEENidS38NATnULpTJ4tQhnXh48pH8dPcpTRZ3SWYa0APgXH8+YXBqfeLHi84FLfjqxuO5cUJ4g7lG9GnHsX4+FCIp1MLuiN5tG0v0c248no+vG8Mn140N8JrWcS2iOMdQw30e2aN1k+cQ2ICaoSF+Y/LFGOvbdG5WBg9OPoJrT+zHF38YF9Y5z/dQfReoaHzz0Tr0AEx94uBiyFrdEjvRrIKIluyMwP6emnOmSghumUBfXrliJFW1Dq59+YegjjuyRxsWbSqPSB4a9HSaMK19Xja/P6l/2OfMiPOxJ2FFJxFpLSJviMgKEVkuIt5nGUpg367ZGessqCT0zbQTmHfziR73hdOAdv6I7gD0KfQ8BfDzlx7d+Pi1K4/h2UsOPnduqyjMzw56Wb/sjPS4WGnpPxcO5/BukZ+7J5wPvkRoFH0I+NgYMxAYCiwPP0vxa2kQ06kq5U/X1i3cVnYKdMrcyUdZQdt19k+As4/sRsn0iR6X1gMY06+w8XH3ti04rn+hx3Q3nDygSfdDgHwvXUBd6xNCrV64/fRDGdDRf+8lf3o2w3TGiSDkgC4irYCxwJMAxpgaY8yeSGUsXlRU1TY+9tRbQilfIr3wyfT/GcITvywKuJdOdkYaFx3jfbUiT6VG5xy75v/c4d0Duq4nOtdM8wunhN4HKAWeFpEfROQJEXFbK0xErhCRYhEpLi0tDeNysXHyA3NinQUVgtOGdErKEZOTj+7B+EMDn5tm5b2nctck/+uJBlodcOOEAdwzaXDA1wc4w+519MrlI72muXpcX47rX8gvhofWLTRaQqk26d0+/CUUAxVOQM8AhgGPGWOOBPYD01wTGWNmGGOKjDFFhYWev9rFs63lVbHOggrBv6YO582rEndBh2hMteqJ6xeKE1xGyuZkpnPhMb38n8f+/dDkI5hm90zytAh4gw6tcnj2kqPjov49ko7u3ZZBncOvQgpUOAF9E7DJGNPQBeQNrACflH7ws2iBUp54quP2JdgamnADv7/jrzqur8dVk1zXkfWW7fycDNLC6BnibdqGRPF0mCONgxVyQDfGbAM2isgAe9OJwE8RyVWccJ4u17XxSqlANFc314bG0yPsVZgmHx163bYvaWlCu7xsnryoiNlOfbh9rSHrLNyeHVkBduu81mXqgXjhvIxhpLqG+hLuX9tvgRdFZDFwBPDn8LMUP5Zsjmy/WKUirXNBDiXTJ4Y8D3vDNwhPgdc5AJ04qGNQdcGRbgz25/cnD2h8PKJ3W/q099xlMxhneBhx/L9276J47Y4eVkA3xvxo148fbow5yxizO1IZiweTHv0m1llQCSycUcWBlmzDLQG7TlAmBF7t09B1kiCOaeApWHpy6eg+PHfJ0fz3N6MCPverVx4TcMneF0+9gw7v1pqS6RNDGuQWaJfUcCR2BVWUXDBSF79VwVn759NC6osd6CGRKgCHU5Ke/ovD6dAqh4c/c58pcWCnfD5fsYP2Xha5vm3iIN5btMXn+T++bkyzj5D1ZViPNrTMSmd/TX3M8hAsDegBuPtM/92+lHIWTkNgMvj9Sf05cVAHhnZvOlpz3s0nsq+6rlmqZL6+6XifywQGKy1NWHb3BHpN+yBi52xuOjGJF869E1L9n1Mlv1Ab7BoWAHHt3piRnsbwnm3d0ncqyPG4IHYkdGuT6/UbQaDGDQisa3UoVV2J0CiatL5eXeY/kVIxFvH5QYIsuwztbtUpR3LelMwM35lozvbWU0Jo97jnrMPI8tGbKRp15w00oHvxq6cXADBJl5tTMeCvNJfo3xl9BeXHpg7nqnF9vc7pMt/LhGbBCrRh1p8LR/Zk3h8jk6dwaUD3ot7ug37nGcENc1YqLHGwXNBlY/rQIT+bEwYGPsVAsHx9s+jeNpebJgz0Ws/eoVVga8n688/zI7dYTU5mfITS+MhFHNMBRSoVXDG2DwCDOrWif8d8vrtlvNfZGlNBpwh8aMy4cHgEchIc7eXigfMIUaXi0Yg+7QA4M0JVgicM7BiVpfUiZUDHfFZur2iWcwfyOgzv2YYNuyo97nvzqmNjNjGcltA9WLZlb6yzoFKcv8bO3u1bUjJ9IqNSZK1MV+9cc3Cg0VlRauf63mk+p7/8z5Am+xrerxaZ6TGd5VMDugcPzPoZIKzFhJUKRexr0BODc3/zBydHZ+H2tk7Vr976u8e6CUSrXDz4fMUOgGbrL6uUSnyXjOrNxMPja9F4Deg+dMiPTGu6UvHmltMGsaBkV6yzkdBuP+PQWGfBjVa5KBWHmrtZ/vKxfZjxy6Jmvkr8+8Wwbjwy5WCVTaJ/K9eA7mL7XmuFoh4hzKamVLhiXQcbLW1yM+nXIY/HY/yh8o/zhnL64QcbVY/q5T5dQaREfFSvB1rl4uKl+RsAUrb3gFLRkJGexszfHxfrbESMz1gdxQ9pLaG7ePyrtQBcMqpXbDOilEo4sf6CpQHdRaU993Gi16UppVKPVrm4KOrZhjSRqC+hpRTA1BE9+WpVmS6qokKiJXQnxhhWbq+gfyctnavYKMzP5s2rjtUusyokWkJ3srW8ioqqOq/TdiqlYuetq49l8+4Dsc5G0K47sR+rtlcwtn9gi2eEI+yALiLpQDGw2RhzevhZip2GOVz6a0BXKu4M69GGYT0OzpNy04SBjO0fH73RcuxFqa89sZ/bvn4d8/n0d9Hp0ROJEvr/AcuB2K3mGiHT3lwMQM92LWOcE6WUP1eN6xvrLDTKSE+Li9kqw6pDF5FuwETgichkJ7Z27q8BrHUPlVIq0YTbKPogcCPg8JZARK4QkWIRKS4tLQ3zckopFV3nFXWjX4J0Yw65ykVETgd2GGMWisg4b+mMMTOAGQBFRUVxu3JEXb31mTSyT/MN/VVKJZ6/njM01lkIWDgl9FHAmSJSArwCnCAiL0QkVzGwrmw/AO3yUnfZLaVUYgs5oBtjbjbGdDPG9AImA58bYy6IWM6ibK0d0C8f0yfGOVFKqdDowCLb6h37AOhbqD1clFKJKSIDi4wxXwBfROJcsbKmdB+dWuWQn5MZ66wopVRIdKSobc2OfTohl1Ip6vPrj2N/dX2ssxE2rXLBmsNlTel+DehKpag+hXkM6VYQ62yETQM6sG1vFfuq6+irAV0plcA0oKMNokqp5KABnYMBXatclFKJTAM6VkBvlZNBoQ4qUkolMA3oWAH9kA55ukqRUiqhaUDH6oOu1S1KqUSX8gF9T2UNZftqNKArpRJeygd0bRBVSiULDegNAb1Ql51TSiW2lB/6v3rHPrIz0ujapkWss6KUioIuBTnsq66LdTaaRcoH9DWl++hTmEd6mvZwUSoVfH3TCbHOQrPRKpfSfTpCVKkUkpYmpCVpAS6lA3pVbT2bdh/QBlGlVFJI6YC+pnQfxmgPF6VUckjpOnTtsqhCYgyIgMNh/Q50hHHlLsi1FyHfuxVatIHMHPd0DgdgoLYSsvOt64Hn6zjqIS394PP6OjD1kJYJaXZ5rWY/SPrBa+3faV07Lc06tzEHH9dVw77tkN8JJA3S7QVfdpdY58jrAI46qN4H+R2t3/U1B+8L4MAeqNoDLTvAzx/BoWcdzGPDa9egtsr6nZlj5b1mH2S2AMS6j90lkNfRykvVHmt7zX6o2Ap9T4DqvbBnI3QcDNUVsOZz6Doc9u+wztf9aOv8VXus43IKICvPei92rbHS7tlg5e/nT6DHSEjPggO7rbQF3WDbUmjbB3LbwYa5ULkT+k+w8vbuNdDvZOt1G3Aq1Fayt6qc/C7DkeyWUL754L0XDvD3FxK2lA7oy7bsBaB3+zirQzfG+kPNbmX989Tsh52rrT/KTQuhfCP88Dx0ORK2/BDxy9cAmUDDv115Whq5Dgfp9ravWuRQK0K3ujra1dXTxuFgfWYG7eqtQORAWJuZSa/aWhblZJNpDB3q68kwhipJo7WjnmVZWexOT6dzXR3/zc/jqANVHFJbSz3QyuFgY0YG+cawLT2dTZkZPN66gP/dW8GKrCxGHzjAo21aU1hXx9DqGg6trqGwvp5ZLXP5MtfqrXTu3gp2p6czq2Vu432dt7eCA5LGjzlZbM7IwBHEVA8tHA4OpKX0F9rQLbyzec77jZftxWGcc/nDgaUrvs36nQmUPG893vSaz0Nmd5lE+5PuDT1vARDT8OkfBUVFRaa4OJxXO7J6TfsAgJLpE5v/YsbAlu+tgPzRDX6Tl6cJq7Ky6Flby+zcXNrW17MxM4Nt6Rn0qq3lsTYF7E5P93sepVR8OCSzNW9P+SqkY0VkoTGmyF+6lC6hZ2ek0bNdrv+EoVj+Pnx+L5Qub9y0IiuTJdlZ3N27R/NcUykVt0IN5sEIOaCLSHfgOaAT4ABmGGMeilTGmlttvYM6h+GkQztG7qSz7oKv72d7ejrje3SFPCAvOYJ3mqThMA4uOvQiahw1fLf1OzrldWJZ2TImD5xMnaOOLzZ+wf7a/Wzdv5ULBk5haIcjaZXVildWvsLFh13M0rKlDGo7iHYt2lFVV8Xa8rUMLRzKnE1z6NyyM4PbDyY7PZtaRy3bK7fTMbcjdY46stKzqKqrouxAGT3ye5Celk6LjBaUVpbSIbcD1fXV7KzaybKdy8jLzOPwwsPJTs8mQzKoddSyoWIDfQv6ArB1/1bSJI2KmgoOaX0IVfVV5GbkIiKUV5dTkF1AraOWLfu2IAi5mbkUZBdQWVtJQbb/JcqMMX5n7aytryUjLcNrOodxYIwhPc39G5gxhqr6KlpkBDYQrs5RR52jjpwMD3X1KumEXOUiIp2BzsaY70UkH1gInGWM+cnbMfFU5bJyWwWnPDiHhyYfwaQjuoZ+ImPgrtZ83DKXGzq0j1j++rfpz7n9z2Vg24HkZebRJqcNLTOtuv70tHQEISMtpb9gKZUymr3KxRizFdhqP64QkeVAV8BrQI8nK7dXANC/YxhzuMx7jBOXPcyOAKpQnj/1eYYWDtU515VSzSYiRTwR6QUcCcyPxPmiYfX2CtIE+oQyStThYPF9HZnapRNkuL+ExRcUk52uqx8ppaIr7IAuInnAm8B1xpi9HvZfAVwB0KNH/NQnr9qxj17tWpKdEWRPkbpqhrxYBF06Ndk8+7zZtG8RuSoXpZQKVlgBXUQysYL5i8aYtzylMcbMAGaAVYcezvUiadWOEFYpqtjGkLdOarLp+wu/JzMtM4I5U0qp0IQ8UkKsyuAngeXGmPsjl6XmV1Vbz7qy/cHVn1ftdQvmSy5aosFcKRU3whn6Ngq4EDhBRH60f06LUL6a1bIte6l3GAZ2DjCgO+oZ8uqoxqePnvAISy5a0ky5U0qp0ITTy+VrDo4OTyjryvYD0Kd9YFUuZz45GLKskvgzpzzN8E5+ew8ppVTUpeTkFCu37SUrI43+Hf0H9Mcf6M46O5g/Ou4hDeZKqbiVkgH9p617Gdgpn4x037e/af1XPNy2NQBXDZjK2J7Ju9KJUirxpVxAN8bww4Y9DOrUynfCmkpO/eLqxqdXj5zWzDlTSqnwpFxAX71jH5U19WRl+L71BW9f1PhYG0CVUokgJQM6wLgBhV7TlN/dmkuqVtAlPZcFUxdEK2tKKRWWlAvoizaVk5kujDrEy6jO/17N6J7dAfjDqLt0ljqlVMJIuYD+48bdHNq5FTmZHob8b1vK77d80vj0pN4TopgzpZQKT0oF9HqHYfGmco7s0cZ9p6OenTPGMNNesuyb872tb6WUUvEppQL6ks3lVNbUc2SP1u47727L39tZgf78gefTKstPLxillIozKbVCwpNfrwNgmGsJ/c4CHmhTwPt5Lbn0sEu5bvh1McidUkqFJ6VK6PuqaslMF7q3dVpH9M4C5uVk81Rra3mxK4deGaPcKaVUeFImoFfV1jN7ZSm/GNbt4MZXL6QOuLyzta7ouf3PDXitRqWUijcpU+Xy+Jy1AAzvaVe3bFsCy9/lSHv5uLuPvZuz+50dq+wppVTYUqaE/uHSbQBMOKwTrPgA/j2aIU5rgWowV0olupQI6Adq6lm+dS9H9WpDvqMCXpnCkwUH50Kfe/7cGOZOKaUiIyWqXF4r3gjAlaO6w/2DmpTMZ54zk7ysIJeiU0qpOJT0AX3zngPc8e4yCvOyGP/mYU2C+fwp88nNzPVxtFJKJY6krnJxOAyjpn8OGL6rO6dJMJ993mwN5kqppJK0JfTXFmzkxjcXk4aDP2f9i8Odgvm8KfNomdkyhrlTSqnIS4qAXlfv4M3vN3HTm03nLR+dtoTx7R/hT+0Ojgxd/MvFiCTkUqhKKeVTwgX0qtp6bnpzMe/8uMVtXx6VXJA+i+szX+H2wna8n9eSRVjB/Jaj/8jkQedHO7tKKRU1YQV0EZkAPASkA08YY6ZHJFdeVNXW88sZcxi1azp/6rCUD7K6siT/QJM0LwIv0qPJthdOe4GhhUObM2tKKRVzIQd0EUkHHgVOAjYBC0TkXWPMT5HKnLOK/RXc+8Id5KZ9yBM9MoG2wAGv6c/scybXH3U9bXPaNkd2lFIq7oRTQj8aWG2MWQsgIq8Ak4CIB/Sbn57ELNZQlSO0r8/kuNzu1OZ34phuY2mR0YIx3cbQuWVnrRtXSqW0cAJ6V2Cj0/NNwAjXRCJyBXAFQI8ePVx3B6Rjy54cU76Po7uO4X/H30hmlnY3VEopV+EEdE/FYeO2wZgZwAyAoqIit/2BuO68h0M5TCmlUko4A4s2Ad2dnncD3LueKKWUiopwAvoCoJ+I9BaRLGAy8G5ksqWUUipYIVe5GGPqROQa4BOsbotPGWOWRSxnSimlghJWP3RjzIfAhxHKi1JKqTAk9eRcSimVSjSgK6VUktCArpRSSUIDulJKJQkxJqSxPqFdTKQUWB/i4e2BsghmJxHoPacGvefUEM499zTGFPpLFNWAHg4RKTbGFMU6H9Gk95wa9J5TQzTuWatclFIqSWhAV0qpJJFIAX1GrDMQA3rPqUHvOTU0+z0nTB26Ukop3xKphK6UUsoHDehKKZUkEiKgi8gEEVkpIqtFZFqs8xMOESkRkSUi8qOIFNvb2orITBFZZf9uY28XEXnYvu/FIjLM6TwX2elXichFsbofT0TkKRHZISJLnbZF7B5FZLj9Gq62j4352oNe7vlOEdlsv9c/ishpTvtutvO/UkROcdru8W/dnqZ6vv1avGpPWR1TItJdRGaLyHIRWSYi/2dvT9r32sc9x8d7bYyJ6x+sqXnXAH2ALGARcGis8xXG/ZQA7V22/RWYZj+eBtxnPz4N+AhrdaiRwHx7e1tgrf27jf24Tazvzel+xgLDgKXNcY/Ad8Ax9jEfAafG6T3fCfzBQ9pD7b/jbKC3/fed7utvHXgNmGw//jdwVRzcc2dgmP04H/jZvrekfa993HNcvNeJUEJvXIzaGFMDNCxGnUwmAc/aj58FznLa/pyxzANai0hn4BRgpjFmlzFmNzATmBDtTHtjjJkD7HLZHJF7tPe1MsbMNdZf/HNO54oZL/fszSTgFWNMtTFmHbAa6+/c49+6XSo9AXjDPt759YsZY8xWY8z39uMKYDnWWsNJ+177uGdvovpeJ0JA97QYta8XMN4Z4FMRWSjWAtoAHY0xW8H6gwE62Nu93XsiviaRuseu9mPX7fHqGrt64amGqgeCv+d2wB5jTJ3L9rghIr2AI4H5pMh77XLPEAfvdSIE9IAWo04go4wxw4BTgd+IyFgfab3dezK9JsHeYyLd+2NAX+AIYCvwD3t7Ut2ziOQBbwLXGWP2+krqYVtC3reHe46L9zoRAnpSLUZtjNli/94BvI311Wu7/fUS+/cOO7m3e0/E1yRS97jJfuy6Pe4YY7YbY+qNMQ7gcaz3GoK/5zKs6okMl+0xJyKZWIHtRWPMW/bmpH6vPd1zvLzXiRDQk2YxahFpKSL5DY+Bk4GlWPfT0LJ/EfCO/fhd4Jd274CRQLn9FfYT4GQRaWN/tTvZ3hbPInKP9r4KERlp1zf+0ulccaUhqNnOxnqvwbrnySKSLSK9gX5YjX8e/9bt+uPZwDn28c6vX8zYr/+TwHJjzP1Ou5L2vfZ2z3HzXseyxTjQH6zW8Z+xWoVviXV+wriPPlit2YuAZQ33glVv9hmwyv7d1t4uwKP2fS8BipzOdQlWA8tq4OJY35vLfb6M9bWzFqskcmkk7xEosv9h1gCPYI94jsN7ft6+p8X2P3Znp/S32PlfiVPPDW9/6/bfznf2a/E6kB0H9zwaqzpgMfCj/XNaMr/XPu45Lt5rHfqvlFJJIhGqXJRSSgVAA7pSSiUJDehKKZUkNKArpVSS0ICulFJJQgO6UkolCQ3oSimVJP4fIEJvp6MGTZ0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "continue_train = True  # set to True to continue to train\n",
    "num_big_steps = 10     # number of small steps\n",
    "num_small_steps = 200  # number of big steps\n",
    "if continue_train:\n",
    "    for k in range(num_big_steps):\n",
    "        for j in range(num_small_steps):\n",
    "            dXY_list = np.append(dXY_list, np.zeros((rep, 1)), axis=1)\n",
    "            dX_list = np.append(dX_list, np.zeros((rep, 1)), axis=1)\n",
    "            dY_list = np.append(dY_list, np.zeros((rep, 1)), axis=1)\n",
    "            for i in range(rep):\n",
    "                minee_list[i].step()\n",
    "                dXY_list[i, -1], dX_list[i, -1], dY_list[i, -1] = minee_list[i].forward()\n",
    "        # To show intermediate works\n",
    "        for i in range(rep):\n",
    "            plt.plot(dXY_list[i, :],label='dXY')\n",
    "            plt.plot(dX_list[i, :],label='dX')\n",
    "            plt.plot(dY_list[i, :],label='dY')\n",
    "            plt.title('Plots of divergence estimates')\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())\n",
    "    display.clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save current results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "overwrite = False  # set to True to overwrite previously stored results\n",
    "if overwrite or not os.path.exists(chkpt_name):\n",
    "    minee_state_list = [minee_list[i].state_dict() for i in range(rep)]\n",
    "    torch.save({\n",
    "        'dXY_list': dXY_list,\n",
    "        'dX_list': dX_list,\n",
    "        'dY_list': dY_list,\n",
    "        'minee_state_list': minee_state_list\n",
    "    }, chkpt_name)\n",
    "    print('Current esults saved.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the ground truth mutual information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground truth is 4.982193620464953 nats.\n"
     ]
    }
   ],
   "source": [
    "mi = g.ground_truth * d\n",
    "print('Ground truth is {} nats.'.format(mi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply moving average to smooth out the mutual information estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "mi_ma_rate = 0.01            # rate of moving average\n",
    "mi_list = (dXY_list-dX_list-dY_list).copy()    # see also the estimate() member function of MINE\n",
    "for i in range(1,dXY_list.shape[1]):\n",
    "    mi_list[:,i] = (1-mi_ma_rate) * mi_list[:,i-1] + mi_ma_rate * mi_list[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the mutual information estimate after different number of iterations. The red dashed line shows the ground truth, and the green dotted line is the number of iterations where 90% of the ground truth is reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VGXawOHfmwJJCIQSOsSAAioIkSYoCkhvsgoiugioC36WVVR2VRZs2FZxFRdXBaUoKKCoiHREUJoUCV06SgkkARJSSXu/P85kMpNMkkmYyUnOPPd15crMqc+cTOaZtx6ltUYIIYTIz8/sAIQQQpRPkiCEEEK4JAlCCCGES5IghBBCuCQJQgghhEuSIIQQQrgkCUIIIYRLkiCEEEK4JAlCCCGESwFmB+AoPDxcR0ZGmh2GEEJUGDt27IjXWtf2xrHLVYKIjIxk+/btZodRYZ1MPAlA47DGJkcihCgrSqk/vHXscpUgxJW5/9v7AVg3ep25gQghPE5rTXaOJsC/7FoGJEFYyMTbJpodghDCS16Yv42tR+JYOWlAmZ1TGqktpGfTnvRs2tPsMIQod/b8cZ6XF27H07NXX0hO56UF20lOz/TocV3ZeiQOgIysbK+fK5ckCAs5dvEYxy4eMzsMIcqd8Z9tYdPBcxw5e6lE+525kEKfyUv5/fRFl+vvffdHNh86x5C3V5XouAs3HaXP5KWMm7nR5fodx+LYdeK8y3Wf/vi7/XHMxdQSnbekJEFYyIOLH+TBxQ+aHYYQVyQh5TJ//3QD8ZfSAfhlfwz7Tl7wyLFXRp8s0fbfbj0OwHNzf7UvG/j6cuauP3RFceR+yB84neBy/YR5W/nn51tcrvtu6wnAaJMYPe2nK4qjOJIgLOTlbi/zcreXzQ5DWExmdk6Znu+t76I5dCaRlxZuJyMrm1cX/cbTszdz6IzrD9OSSM8sWD3z6PRf6DN5qcvtz1wwvqGnZWTTZ/JSsrJzyMzO4fOfD19xLLk+Xr2fJ2du5I+4JIb/Zw1xl9Lc2m/N7tMei6EwkiAspGtkV7pGdjU7DFFCRu+Usv0QdteqXScZ+PpyzpagKiNHaxZsPML5pPRSnTPQ1kundrUgft4fY1/+90838tqi34rcV2vNjmNx/HIgb7+xH623P64XFlxgn6PnjGqnS6kZBdZVCnD+iHSMx5UNB2JYvetUkdvk982W4/x+OoHvtp7gYsplp9LK2YTCr/uU73eV6DylIQnCQg7GH+Rg/EGzwxDAzB9/Z+T7a93atu+ry+j/2nLOJqSyZvcp0jOyPBpLbGIa7/6wm9hE52+ml9Iy+GZL0W1W73y/G4BtR+M46MY3+F8OxPDZukPMXHuQ5+f96rTu/z7+mT6Tl/L24ugiG1q3HI4FoFKAP28vdv4Q/Hl/DH0mLy20ymneL0eYMG8rr379G30mL+Xo2UT+iEu2r88oojT0yPRfnJ4fiUlk08FzTsv+/V20/bHjN/0B7SK4kJzO5K9/Y8r3u4psG1i+80+Xy5f9Ziw/dT7FvmzUf11XIXmqyq04kiAs5OEfHubhHx42OwwBLNh0lHOJxVcV/OLwjfS5ub/y9uJdvLRwh8ttv9p0lDv/vbLEPXHuf38tK3ae5P58CeuJTzfy8eoD9uqV1MtZaK0Z8+F6Vu86xbhZeQ2o05bv5YlPNzLhi61FnuvVr3/jyw1HAPgjLplnbfXoUxbv4nhsEmBUjRyOSSw27nX7zhS67unZm10u/zxf28CjMzY4PV+w8SgHTuU1OGc5JIz4fCWeaSv2FhnfiKnO1/N+h+cvLSh8wO97P+wp8rj5vb04ukAVWPRx1w3YnibjICzk9R6vmx2CT0lMzaBacCBKKaflSWmZTtuEhVRyuf/J+GRedagyyf3WufN4PL8di6dt03Cn7T+xNWxmZOVwOTMbfz9FlaBAYi6mMnPt7zzR/wa+23qcGqGVGdjuKi4kpxMaFFho/I7fcudvOMKsnw7SsGYVTl9IKbT6YsfROLTWBV4zGFVL+UWfOO+yfj87p+C2sYlpBZJYUbYcOken5nUBmLp0D+FVg9zaL+5SOtfZHo9xqH4C6DN5KSsnDUBrzYFT7rd5ZGTmkOXwmk7EJfHFL4eZs+4QD97egntuucbtY+Xnqq1h/sYjpT5eSUiCsJCbG99sdgg+Iz0ji2HvrKZTszq8PLyDffn4OZvZ82de8X/YO6sLHdj0tw/Xu1wO2KtnPn74NiLrVHUqNew6cZ5J87cBsHLSAHtPFsf68QOnLrr8YMnO0ayMPsntrRowoF0ES3cY1RqzfjKqJk9fSCmwT35bj8RyU7O6BZZPLKZ04SgtXzXaiKk/EnepZG0WLy7YzspJA7icmW2vnnHHa4t+o0rljrS7ura9EdqR1ppLaSUb17B6t3O7w3WNqjNnnVGambn2oFOCqBsW7FbpsigZWWXTZiVVTBayN3Yve2OLLhaLKzP358P0mbzUXtWSW18OxoeeY3LItf1oHP1fW0Zmdg7nElLZnK9euygPf/wzOVo79STKTQ4A/1ni+pt+YT1c1uw+xdSlexj875VuJQNXdhyNd738mOvlrrwwfzt7/sirJnEnOeRvMM41/D9r3D5vriXb/2D3H66raTYfOsewd1bbn19VO5SVkwZQo0plt48fXKnw796tI2u5H6jJpARhIY8vexyQuZi8JfVylr2Oe9/JvHrslMuZVKkcWOgHzr9syWTg68sJDQogOb1kjdAvLtjOM4Nau1y3MrpkPWb+s2S3/XFp67EXbzvB4m0nipzyISykEokuegU5Gv/ZFlZOGuDWKOTcc+Wvrhr137WkutGoP+K2ZhyJSbQn9M2HzrHtSKzLbV/O1wb04dhbAahVtTIXUy4Xey6A31wky+wco3dXSXs5Obrl2nps/P1sqfcvKSlBWMjbvd7m7V5vmx1GhZecnsnbi6NJcfjg2nfyAne+tdLl9ne9ZYyiDfDL+3d6+Z72hRy74IfZion9i4xn6+FYe3VFaQ27+Wqn5+HV3Kuvd9SwZhWn519tNkYDf+YitoXP9HJ6HlzJnyXP98UvX9tFn8lLXY5Cdje+swkFq2puu74+r9/X0f781Xs7MOK2Zky6u53Tdlku2kFc8bf9XV8cZvxNB7W/yr5u1mPd3DpGn8lLeeLTDSX6O37/XN8Cy/45uI3T87s7N3X7eKXh1QShlDqhlNqjlIpWSsk83l7WoWEHOjTsUPyGokgLNx5lze7T3PX2KtbsPsXHq/cX2msm19z1h+zfhP/SMZKbmtVx+3z5G3xru/hwLEkde34Lnu7JkE5NnJbFu1nf/9Kw9nzzj968em8HZub7MPxkjdFoPu8X50Fjcx7vDjgniS+f6kmlAH+WF5MMc43peV3xGxXiX0Pa0u7q2kwa2paHe19Ph2vqoJTCz69gwzrA9Y1qcN+trhuR5/y9u/1xnbBgVk4awAPdW9iXNahZhacLKd3ll3+aj6/H93a53fg72jD78e5UDvSnrsO4jWf/EkVQvqqrv13BdXJHWVQxdddau185KUot+qzRRzuqXpTJkVRMF5MvczkrmwWbjtqX5e+HX5jPfz7MVbVDAbi789UopXjx7na8/JXrLqu5boioCcDj/Voybfk+41hP3M6Rs5d4/JMNRe7n2N7x4dhbqVc9hI2/n3XqgXTrdfWoXkTdeeUAPy7bGjw/fvg2ABrUDOHQmUQiwkOpZuuB1eGaohOeYw+mejVCAKOaKbdBtqg6eVeqBgcy+/HuHDqTQC2H3kmP9m3Jr4fO8Wd8coF2i6cHtaZPVN69ULpcV99pvZ9SvHpvByZ+abThNKpVhVPnU/jP6M6A0fjvWHUIUK96SIHYQioHcFXtUEZ3MxJFn6jGTlV3JXmNIZUCnKrIfpjQzz5QEGD237vT79VlALQxoe1C2iAsZNyKcYA12yBy656/GNfD6QOjNDKzc/ho5T7uu7WZ07GGv+teY2f9GiG8M6ozVYICeX/pHn7cYzQI5w7ICg02upbefG29Yo/11shOAPS9MYJpy/fx7F+iUErRrH4YKycNcNlF1E8p3hhxEwNfX25fVrd6MCGVA+jVphFnE1KZa5sKYuLQvGqVYTdfzUKH5AfQODyUN0bcRFJaplMVUitb4nLXGVuD9103OZdUPnvidreP8bce19IyoiYfr9rPjU3C8VOK+jWcP6AHd4hkcIdIzlxI4YEP1jmtc0wOhelwTR16tm7Imt2nOXU+hZtb1LWX4G5qVscpQbS/2vVN2pRSTP8/5xkLhnRqwqItx4msXZUTcUnuvFwAZjzSle1HY9l48ByDO0Q6JQcw/tY/TOjHpdSMK37fl4a3E4QGVimlNPCx1np6/g2UUmOBsQCtK1eGbt2cNxg2DB59FFJTob+L4uno0cZPfDwMHVpw/SOPwD33wMmTcP/9Bdc/8wwMGgQHD8LDLgaZTZwIPXtCdDSMG1dw/euvw803w6ZNMGFCwfXvvQdRUbBmDbz6asH1H38MLVrAkiXwzjsF13/+OTRuDAsWwIcfFlz/9dcQHg6zZ/Pe17YRo7O75a1ftgxCQuB//4OFCwvuv26d8XvKFPjhB+d1wcGw3PYhNHky/Pij8/patWDRIuPx88/D5nzVMI0awdy5xuNx44xr6Kh5c5hue0uMHQuH8tXPRkXxzfAn+NWhp9B97/3I4p/eIygnk5zOnflu8Bju6tQUhgyB80ajqwb69vgHPTjPPyeNNHbs1w/SjPrq2dd05YerOvLDjpJX23SJPcgz65YT8nUmjB7NyMHD7AkiV1CvHsaDRx4BQgsc45t1U1lb73paDu5h1MkfPEjgww+zEsDxEk+cyIKne3KPQy+dORunUy89Eda8xQfPv8JjW40Poyp98qpz7gfm3T6evlfZPvBt772FPf5RIJaWjWtSbc1Kqrn53vtLs+58F1GwfeWh/xlddhO/+wFef9RYWMh7z7/702T7+Tvt/+62eVz/ozEwbirAV0W/94LqNoJW99qfDz65A7o5tL8V8d77pds48DeSeHClAPt7r3mNCGh7j327Z++0lcRHjIBT+RqWO3eGN94wHg8ZwrX+teCGwQQcPUJUVjrRNY12iqgmtVx2BujvbySi8GpB9H1qJH0BPnLYwOFzL7B/fxzLDn+L6MAnzbrxzT97G597XuTtRupbtNZtgX7AY0qp2/JvoLWerrVur7VuHxhY+KAeUbyo5FCikgt+IFVUj4R14uPVB4jON+3x4O7jSA6oTL+QLny8+gBvfRfN3ip1OVzVqAZJDjCqVH6kln3aiuHtH+CTa7pyuGodvr6qI6U1ac/3hGTnNV7XCHWuvqmb5jxCeMrITgT4Kapm5jWmVsnOYNDpaJoWnBaogOpVKrNyUEP783rpeccPtH3zDQssWLe+dO07PNm6utOyay4V7F4bFOhfYFlRHjn8Eyt/LLwjRKf44gdwPXLIeTBcJ78krr9U+KhpV6plOVcvLW7crpAtC+pzJm8ks2Nyj0jJe5/1Tz9FtWDXAxxdUbYattigqvbk0KlZHcb2vL7Ati0TTvFkpdL3RLr7z22sbBxHlcre/7xUnr6BRqEnUuolIFlrPaWwbdq3b6/lntSlt+20UbdqlYbqwmbYLMqKif2ZseYAi7Ycty/r1rJBkdM2FOZ/Y27l3R92M7JrcybN38ZXz/Sy18k7Ono20WlKB1fdPx1fS2nuCJa7v+O+J2KTePjjn2lUqwqfPtqt2GMkp2dy37tr+Pf9nRg3axMAD/e6ziiBlVBGVjbpGdks3HSUrzbnzec0qltz7ru1WbH7547chtJdDzAGtPW11c83qVOVjx4u8P3TpT/ikhj70c/257nndzzePwe3oUfrRm7H8tZ30fZkU71KJRJSMri9VQP+r09LpzEVYPSqKq5NpySUUju01q67zV0hr5UglFJVlFJVcx8DvQEZxeVF/1j9D/6xumA1QnmWlpHFLtt0DI4jgU+dTy5ir8L1fXWZU3KAouf0cfTuA3kj0Wc+1o2r61Vj2t+60LFZHVZOGuAyOQA0rVut2GNPffAWAOY/Vbo7/n329+72RuRcub2dHrz9WreOERoUyPfP96OOQ8+Y/m0jShVPpQB/qoVUolvLBk7Lb7u+fiF7ODtZyr+vI6UUr9/Xkeb1w/hgTBe397uqdlX748EdIp2Ol6tdIe0PhansUBIbaut6enW9MEIqO9fiz3qsm0eTg7d5sw2iLvCt7aIHAF9orVd48Xw+b1r/aWaHYLcy+iT/WbKbRjWr8GkhfcUTUzOcvl29tug3XlsEXz7Vw16n7aj1VTXZ/YfnZ7H85+A2RNapytX1wuy9bvL3+S+KUoqXhrXnpYXbnfrIO7q2YfUrupdwXRe9aaoEBZbqmLWqBtH1+voM7dy0QLfJkso/JUWjWu5VceZu91AP95JbYdpdXbvEH+aO8s9wm6uonl+u9LuxMct++5Nn7mhN91YN0RoGdzQane+/rZn9/hENSvC+Kg+8liC01seANsVuKDymVZ1WZXq+lMuZ7DgaT6fmdagU4FyXndvt79SFFJLSMqkaHEi/V5fZu0N+88/ehU71cO+7eQ2S/dtG8OSAG5zWD3tndbGjdIvz/J038sa3O+1zHeWa9Xg3MrNLXu3auUXdMr2Z/JWaMKStR46Tv5eRu3J7O+Xv9VTWerdxrkaa+uAtHI8t2W1JAZo3cP4C4DgwcUTX5vx6JNaptFJRlFkbhDukDeLKbDpp1CuX1aR9RdWrO657uNd1fLz6gNP6qsGBTrOeFmbmY91cfpvP3z4x/6meBbqpVgsOdDnpWkX6IC/vzielc997RkJv2bgG/xldMSaMzMrO4fSFFKfqporKm20QMg7CQib8aHSzLYtxEK/kGwD26tc7mDCkLY98/AuDO0Y6rcufHAC3kgMUnN4h19IJ/cjRmgB/P7JzNIH+fnRuXpfNh4yeOjMf7Ua9GsH0f83oKtnvxsYs33mSJc8XnL5AlF5NWy+uwR0iebRvS5OjcV+Av58lkoO3SQnCQnLvJtcivEUxW5bego1HCarkz/9W7Ct22790jLTfYD3X8on97SNDHfW4oSFPDWrNm9/sZINtMrKFz/Qq9F4KrlzOzOaz9YcY3CHS3hB78EwC2Tma5vXDSErLLNAtVYiKTkoQwi3eTAxgzEY5c+3vTsuKGjna78aIAgki/2RtwZX8WfhMLwL8/fBTiufuupHVu07RJ6ox/oXMnVOYyoH+BebwadEgbyyAJAchSkYShIWsP2H0/Oka2bWYLUtnb757HTw54Ab6t42w3zkrv8g6VfnvQ7fw90+NW1fe28WYEG3lpAEcPJPA7hPnuTvfLKOB/n6l7nophPAsme7bQl5c9yIvrnvRo8d0vDXkP233F87V1dbn3dXAqNwZLps3qM6MR7ripxR9HebKadGgeoHkIIQoX6QEYSEzB8+84mPs/uM8p86ncDkzm49W7S9y2yoO9zv+7O/dGflf49aX+XsJRYSHuj3NsxCi/JAEYSFNa5Tu5iE7jsYx4Yut9Ilq5NYdyqoGB/JoH+ceK3Wrh7Dk+b4lbjcQQpRfkiAsZM0xYxxAz6buT+ew588L9vsru5McHuvbkjsKGfCTf7CcEKJikwRhIa/+bEwn7pggLmdmc+RsIi0bG/P77/3zAtk5mrSMLF5c4F6X4tfv68iEL7by+n0dr2haAyFExSIJwkI+v/PzAsve+GYnmw+dY8rITny54Qg7XNxMHShQvdS8QRhv3d/JficwGX0shO+RBGEhjcPyegnlDoDMHVk8/rMtLvcBqFW1Mk8PamNPEJIMhBAgCcJSVhwxJsvte01fXlyw3elObIUJDQrkc9ttIZc835dyNLBeCGEySRAWobXm+VWv4O+nOPNHZKHJYcqozvyw/Q9+ORDDF+N6OE1rLI3MQghHkiAsYtS0nwi7aNwLeEFc3o3p/3prM0Z2a87kr3YQER7KDRE1uSGiJs9zo1mhCiEqCEkQFnEuIY3KqobTsm//2cd+R6tJd7t/z14hhABJEBVW7hQYOVoz9O1VAMTqXwGoo27iq2d6FbjdoRBClIR8glQw6ZnZJKZctk9r4Sg5dDkNa1Zh5ehXTIhMCGE1kiAqkAvJ6U6348xv/l0LuaqO3ARFCOEZkiAqkNHT1rlcPvyWqxnZrYXMgySE8ChJEOVcVnYOvxyI4c1vo52WTxzSlospl4kIDyWqSTgA3xz4BoC7rrurzOMUQliPJIhy7Ni5Szwy/ZcCywsb6fz+r+8DkiCEEJ4hCaIcc5Ucnv1LVKHbLx6+2JvhCCF8jCSIciZHa77adJSZaw86LR9/Rxtual6HasGVCt03LCjM2+EJIXyIJIhypt+rywosc3fyvAV7FwBwT6t7PBqTEMI3SYIoRw7HJDo9f6J/K267voHb+3+4/UNAEoQQwjMkQZQTGVnZPP7JBgDq1whh9uPdS3yMZX8tWPoQQojSkgRRDuTvrTTjka6lOk5IYIinQhJCCEkQZuszeanT86UT+hHg71eqY83dPReAEa1HXHFcQgghCcIkWmv65muQnvVYt1InB4BPfvsEkAQhhPAMrycIpZQ/sB04rbUe6O3zVRSOyaFSgB/fPdsHf7/SJweA1fevvtKwhBDCrixKEE8CB4BqZXCuCuGxGXntDc/+JYrbb2jokeMG+gd65DhCCAFwZV9Zi6GUagQMAD7x5nkqks0Hz3Hk7CX7c08lB4DZ0bOZHT3bY8cTQvg2ryYI4D3gn0COl89TIeRozUsLt9ufL36ur0ePLwlCCOFJXqtiUkoNBGK11juUUt2K2G4sMBYgIiLCW+GUCxO/2Gp/vHxif/yUZ6fnXjd6nUePJ4Twbd4sQdwC3KGUOgHMB25XSs3Nv5HWerrWur3Wun3t2rW9GI65vv31ODuOxQMw87FuHk8OQgjhaV5LEFrr57XWjbTWkcBwYK3W2if7X247EstHq/YD0PqqmjSsWcUr55mxYwYzdszwyrGFEL7H220QPm/tntNM/HIbAP5+irdHdvbauRbsW8CCfQu8dnwhhG8pk4FyWut1wLqyOFd5cjkzm39/Z9wJ7tqG1Zn64C1ePd+akWu8enwhhG+REoQX3fHmCvtjbycHIYTwNEkQXnI2IdX+2N37OVyp/237H//b9r8yOZcQwvokQXhBWkYWo/77EwBvj+xUZuddcmgJSw4tKbPzCSGsTSbr87Dk9EyGvL3K/rz1VbXK7NzL/7q8zM4lhLA+KUF42Fu2RmkwBsMJIURFJSUID9pxLI5fD8dSOdCf757tU+aD4aZumQrAk52eLNPzCiGsSUoQHjRhnjGVxiePdDVlpPSPx3/kx+M/lvl5hRDWJCUID7mYfNn+uE5YsCkxfH/v96acVwhhTVKC8JDcWVo/GnuryZEIIYRnSILwgIysbH4/nQBAZJ2qpsUxZdMUpmyaYtr5hRDWIlVMHjBn3SEAJg5tizJxltbNpzabdm4hhPVIgrhCKemZfL35GABdrq1naiyLhi0y9fxCCGuRKqYrtH5/DAD/fegWU0sPQgjhaZIgroDWmjW7T1G/RgjN6oeZHQ5vbniTNze8aXYYQgiLkCqmK7D9aBz7Tl7k8X4ty0XpIfpsdPEbCSGEmyRBlJLWmnk/H6ZuWDB9bywf99KeP3S+2SEIISxEqphKae2e0xw4ncDdN19NoL9cRiGE9cgnWynkaM1bi3cB0CeqkcnR5Jm8fjKT1082OwwhhEUUmyCUUs2VUj8qpfbanrdWSk30fmjl1/Lf/gSM5FApwN/kaPIcPH+Qg+cPmh2GEMIi3GmDmAH8A/gYQGu9Wyn1BfCqNwMrz95ftheAJwe0NjkSZ3Pvmmt2CEIIC3GniilEa70137IsbwRTEZyITQLg7s5N8fczv+eSEEJ4izsJIl4pdTWgAZRSQ4EYr0ZVjq3dexo/pRjauanZoRTwwk8v8MJPL5gdhhDCItypYnoMmA5cq5Q6DRwH/urVqMqpHK1Zt/cMbZuGU71KZbPDKeDkpZNmhyCEsBB3EoTWWvdUSlUB/LTWSUqpJt4OrDzaf/Ii5xLTGNWtudmhuDRr8CyzQxBCWIg7VUyLALTWKVrrJNuyr70XUvn1097TVA7wo3MLcyflE0KIslBoCUIpdS3QEghTSt3lsKoaEOTtwMqbrOwcft4fQ+cW9QipXD4HoD+/5nkA3uj5hsmRCCGsoKhPuhbAQKA6MMhheRIwxptBlUc7jsVxKS2T7q0amB1Koc6nnTc7BCGEhRSaILTWi4HFSqnOWmufvxPN2j1nqBocSLura5sdSqGmD5pudghCCAtxp65kp1LqMYzqJnvVktb6Qa9FVc6kZWSx+dA5etzQUOZdEkL4DHc+7T4H6gF9gPVAI4xqJp+x+eA5Lmdmc3s5rl4CGL9qPONXjTc7DCGERbiTIK7RWk8CUrTWc4ABwA3eDat8Wb/vDLWrBdEyoqbZoRQpLTONtMw0s8MQQliEO1VMmbbfCUqpVsBZILK4nZRSQcDPQGXbeb7WWr9YyjhNk5GVzc4T5+ndphF+5eCmQEX5YMAHZocghLAQdxLEdKVUDWAS8D0QCrgzn8Nl4HatdbJSKhDYoJRarrXeUvpwy96ePy9wOTObdk3Lb+O0EEJ4Q7EJQmv9ie3hesDtCYi01hpItj0NtP3okgZotq2HY6kU4MeNTcPNDqVY41aMA+C9vu+ZHIkQwgqKTRBKqerASIxqJfv2Wusn3NjXH9gBXAN8oLX+tdSRmmTbkTjaRNYiKLD83PdBCCHKgjtVTMuALcAeIKckB9daZwNRtiTzrVKqldZ6r+M2SqmxwFiAiIjycW/nXKfPp3D6Qgp/6RhpdihukZKDEMKT3EkQQVrrp6/kJFrrBKXUOqAvsDffuukYs8XSvn37clUFtfVILAAdrqljciRCCFH23BoHoZQao5Sqr5SqmftT3E5Kqdq2kgNKqWCgJ/D7FcZbprYdiaVxrSrUrxFidihueWzpYzy29DGzwxBCWIQ7JYgM4G3gX+Q1MmuKb7CuD8yxtUP4AQu11j+UNtCylp6Rxe4/LjCow1Vmh+K24MBgs0MQQliIOwniaYzBcvElObDWejdwY6miKgd2Hj9PZnYOHStQ9dKU3lPMDkEIYSHuVDFB8yeWAAAfQElEQVTtA1K9HUh5s+1oLMGV/GlVzkdPCyGEt7hTgsgGopVSP2EMfgPc6+Zake0+cZ4bImpWqMn5xi4ZC8isrkIIz3AnQXxn+/EZCSmXOXk+hV5tGpkdSonUCq5ldghCCAtxZyT1nLIIpDzZf+oiAC0bV6zqJbmTnBDCk4q65ehCrfUwpdQeXEyRobVu7dXITLT/5EUC/BTNG4SZHYoQQpimqBLEk7bfA8sikPJk/6mLNKsfRqWAijW9xgOLHwBg1uBZJkcihLCCQltgtdYxtoePaq3/cPwBHi2b8MpeRlY2h84kcn3jGmaHUmKNqzWmcbXGZochhLAIdxqpewHP5lvWz8UySzhy9hKZ2Tlc36jiJYhXur9idghCCAspqg3iEYySwtVKqd0Oq6oCG70dmFn2nzQaqCtiCUIIITypqBLEF8By4A3gOYflSVrrC16NykSHYxKpXS2ImqFBZodSYiO+GQHA3LvmmhyJEMIKCk0QWutEIFEpNRE4q7W+rJTqBrRWSn2mtU4oqyDL0uGYRJrVr5i9l1rUamF2CEIIC3GnDWIR0F4pdQ3wKcZtR78A+nszMDOkpGdy+kIKPVs3NDuUUpnUdZLZIQghLMSdeSRytNZZwF3Ae1rrpzBmarWcI2cvAVTYEoQQQniSOwkiUyl1L8ZtR3On6w70XkjmORyTCFTcBDH86+EM/3q42WEIISzCnSqmB4D/A17TWh9XSjUBLNkKevRsIuFVg6hepbLZoZRKVL0os0MQQliIO3Mx7VdKPQtE2J4fB970dmBmOB6bRNO6Vc0Oo9Se6/Jc8RsJIYSbiq1iUkoNAqKBFbbnUUqp770dWFnLys7hZHwykXWqmR2KEEKUC+60QbwEdAQSALTW0UATL8ZkijMXU8nK0USEh5odSqkNWTiEIQuHmB2GEMIi3GmDyNJaJyqlHJcVmN21ovszLgmAq2pX3ATRuVFns0MQQliIOwlir1LqPsBfKdUMeALY5N2wyt6f8ckANK7AJYjxN483OwQhhIW4U8X0d6Alxu1GvwASgXHeDMoMJ+OTqRMWTHAld3KmEEJYnzu9mFKBf9l+LOvP+OQKXXoAuOPLOwD4/l7L9SEQQphAvi4DOVpzMj6ZG66q2Pd07tGkh9khCCEsRBIEEJuQxuWsnArdgwngyU5PFr+REEK4yZ02CMvLbaCu6AlCCCE8qagbBv2XIrqzaq2f8EpEJvgj3ujiWtETRL95/QBY/tflJkcihLCCoqqYtpdZFCY7GZ9M9SqVqBZSyexQrsig5oPMDkEIYSFF3TBoTlkGYqY/45MrfOkB4NEOj5odghDCQoqqYiqyr6TW+g7Ph1P2tNb8GZdM91YNzA5FCCHKlaKqmDoDJ4EvgV8BVcS2FdbFlMukXM6yRAmi52c9AVgzco3JkQghrKCoBFEP6AXcC9wHLAW+1FrvK4vAykrMxVQAGtSsYnIkV+6elveYHYIQwkKKaoPIxpjie4VSqjJGolinlHpFa/3f4g6slGoMfIaRaHKA6VrrqZ4J23PO2hJEveohJkdy5ca0G2N2CEIICylyoJwtMQzASA6RwPvAN24eOwt4Rmv9m1KqKrBDKbVaa73/CuL1uJiENBRQt3qw2aEIIUS5UlQj9RygFbAceFlrvbckB9ZaxwAxtsdJSqkDQEOgXCWIsxdTqVUtiEoB/maHcsW6ze4GwLrR60yNQwhhDUWVIO4HUoDmwBMO94NQgNZau33rNaVUJHAjRmN3/nVjgbEAERER7h7SY2ISUqlvgeolgNFRo80OQQhhIUW1QXhkGg6lVCiwCBintb7k4jzTgekA7du3L/MbEZ29mMqNTcPL+rReIQlCCOFJXp2LSSkViJEc5mmt3W27KDMZWdnEJ6VbpgSRmZ1JZnam2WEIISzCa7O5KqNO6lPggNb6P946z5U4m5AGQP0a1kgQvT7vBUgbhBDCM7w53fctGO0Ye5RS0bZlE7TWy7x4zhKJuZgCQD2LJIi/tf2b2SEIISzEawlCa72Bcj76Om8MhDW6uI5oPcLsEIQQFuLT94OISUijcoAfNapUNjsUj0jNTCU1M9XsMIQQFuHTd5SLuZhKvRohOHThrdD6z+sPSBuEEMIzfDpBnL1onTEQAI+0f8TsEIQQFuLTCSI2MY02kbXMDsNj7mklk/UJITzHZ9sgUtIzSc3Iona1ILND8ZjE9EQS0xPNDkMIYRE+W4KIu5QOQLiFEsTg+YMBaYMQQniGzyaI+CQjQdSuZo0urgBP3PSE2SEIISzEdxPEJWMUtZVKEHddd5fZIQghLMRn2yDibVVMtapaJ0HEp8YTnxpvdhhCCIvw2RJE3KV0alSpTKC/dXLk0IVDAWmDEEJ4hu8miKR0S1UvATzT+RmzQxBCWIjPJoj4S2k0rFnF7DA8alCLQWaHIISwEOvUr5RQ/CXrlSDOJp/lbPJZs8MQQliET5YgUi9nkXI5i/Cq1uniCjD86+GAtEEIITzDJxNEbhdXK42iBniuy3NmhyCEsBCfTBBxSdYbRQ3Q95q+ZocghLAQn2yDyB0DYaVR1AAnE09yMvGk2WEIISzCJ0sQeYPkrHGjoFz3f3s/IG0QQgjP8M0EkZROWEglKgX4mx2KR028baLZIQghLMQnE0TcpTTLNVAD9Gza0+wQhBAW4rNtEOEWa38AOHbxGMcuHjM7DCGERfhoCSKdVhE1zQ7D4x5c/CAgbRBCCM/wuQSRnpFFcnom4RaaxTXXy91eNjsEIYSF+FyCiLfoGAiArpFdzQ5BCGEhPtcGYcVbjeY6GH+Qg/EHzQ5DCGERvleCsOggOYCHf3gYkDYIIYRn+FyCiMu91agF2yBe7/G62SEIISzE5xJEfFI61YIDqRxorUFyADc3vtnsEIQQFuJzbRBWHQMBsDd2L3tj95odhhDCInyvBGHBGwXlenzZ44C0QQghPMNrCUIpNRMYCMRqrVt56zwlFXcpjWsbVTc7DK94u9fbZocghLAQb5YgZgPTgM+8eI4SuZyZzaU0aw6SA+jQsIPZIQghLMRrbRBa65+BC946fmnkDpKzYhdXgOiz0USfjTY7DCGERfhUG0S8hQfJAYxbMQ6QNgghhGeYniCUUmOBsQARERFePVe8hcdAALzX9z2zQxBCWIjpCUJrPR2YDtC+fXvtzXPF2UdRWzNBRNWLMjsEIYSFmJ4gylJ8UjqhQYEEVbLmy952ehsgjdXCMzIzMzl16hTp6elmhyKAoKAgGjVqRGBgYJmd05vdXL8EugHhSqlTwIta60+9dT53xF1Kt2zpAeAfq/8BSBuE8IxTp05RtWpVIiMjUUqZHY5P01pz/vx5Tp06RZMmTcrsvF5LEFrre7117NKKv5Rm2QZqgGn9p5kdgrCQ9PR0SQ7lhFKKWrVqERcXV6bntWZdSyHik9JpVj/M7DC8plWdcjMeUViEJIfyw4y/hc/MxZSRlU1CSoZlx0AAbDq5iU0nN5kdhhCWsG7dOgYOHFhgeXR0NMuWLSvVMV9/PW/G5RMnTtCqVfn+UuczCeJ80mXAumMgACb8OIEJP04wOwwhykxWVlaZn7OoBFFcPI4JoiLwmSom+xgICyeIjwd+bHYIQnjM5MmTmTdvHo0bNyY8PJx27doxfvx4unXrxs0338zGjRu54447GDp0KA8++CBxcXHUrl2bWbNmERERwejRoxk4cCBDhw4FIDQ0lOTkZNatW8dLL71EeHg4e/fupV27dsydOxelFCtWrGDcuHGEh4fTtm3bAjFlZGTwwgsvkJaWxoYNG3j++ec5cOAAZ86c4cSJE4SHh9O7d2+2b9/OtGlGm+DAgQMZP348K1asIC0tjaioKFq2bMlrr71GdnY2Y8aMYdOmTTRs2JDFixcTHFx+ajl8JkHYx0BYdJAcQIvwFmaHIKysW7eCy4YNg0cfhdRU6N+/4PrRo42f+HiwfVDbrVtX6Km2b9/OokWL2LlzJ1lZWbRt25Z27drZ1yckJLB+/XoABg0axMiRIxk1ahQzZ87kiSee4LvvvivypezcuZN9+/bRoEEDbrnlFjZu3Ej79u0ZM2YMa9eu5ZprruGee+4psF+lSpV45ZVXnBLASy+9xI4dO9iwYQPBwcHMnj3b5TnffPNNpk2bRnS0MR3OiRMnOHz4MF9++SUzZsxg2LBhLFq0iBEjRhQZe1nymSqmvHtRl5/s7GnrT6xn/Yn1ZochxBXbsGEDgwcPJjg4mKpVqzJo0CCn9Y4f3ps3b+a+++4D4P7772fDhg3FHr9jx440atQIPz8/oqKiOHHiBL///jtNmjShWbNmKKVK9EF9xx13lOqbf5MmTYiKMga4tmvXjhMnTpT4GN7kQyWINEKDAgipbN2X/OK6FwEZByG8pIhv/ISEFL0+PLzo9floXfSkClWqVCl0XW5vn4CAAHJycuzHy8jIsG9TuXJl+2N/f39720Fpewo5xuN4XqDIgYb540hLSyvV+b3Fd0oQiWmW7sEEMHPwTGYOnml2GEJcsS5durBkyRLS09NJTk5m6dKlhW578803M3/+fADmzZtHly5dAIiMjGTHjh0ALF68mMzMzCLPee2113L8+HGOHj0KwJdffulyu6pVq5KUlFTocSIjI4mOjiYnJ4eTJ0+ydetW+7rAwMBi4yhPfCZBxF5Kp3aYtRNE0xpNaVqjqdlhCHHFOnTowB133EGbNm246667aN++PWFhrscwvf/++8yaNYvWrVvz+eefM3XqVADGjBnD+vXr6dixI7/++muRpQ4wprKYPn06AwYMoEuXLlx11VUut+vevTv79+8nKiqKBQsWFFh/yy230KRJE2644QbGjx/v1Ng9duxYWrduzV//+ld3L4WpVHFFubLUvn17vX37dq8ce+iUVXS9vj5/73+DV45fHqw5tgaAnk17mhyJsIIDBw5w3XXXmXb+5ORkQkNDSU1N5bbbbmP69Okuexb5Eld/E6XUDq11e2+cz7oV8g7SM7JISsukjsVLEK/+/CogCUJYw9ixY9m/fz/p6emMGjXK55ODGXwiQcResvad5HJ9fufnZocghMd88cUXZofg83wiQcQlGj0DrN4G0TissdkhCCEsxCcaqWNto6jrWHgUNcCKIytYcWSF2WEIISzCR0oQ6SigloVHUQO8ueFNAPpe09fkSIQQVuATCSL2Uhq1qgYR4G/tAtP8ofPNDkEIYSHW/sS0ibuURu0wa5ceAOqF1qNeaD2zwxDCI6ZOnUqrVq1o2bIl7733nn35hQsX6NWrF82aNaNXr15cvHgRgEWLFtGyZUtuvfVWzp8/D8DRo0cZPnx4mcbtiWm8Q0NDPRTNlfGJBHEuIY06Fu/BBLDk4BKWHFxidhhCXLG9e/cyY8YMtm7dyq5du/jhhx84fPgwYEx616NHDw4fPkyPHj14802javWdd95hy5YtjBw50t4DauLEiUyePNmtc2ZnZ3vnxVRglk8Q2Tk5xCamUb9GiNmheN07m9/hnc3vmB2GEFfswIEDdOrUiZCQEAICAujatSvffvstYEybMWrUKABGjRpln7nVz8+Py5cvk5qaSmBgIL/88gv169enWbNmhZ4nNDSUF154gZtuuonNmzezY8cOunbtSrt27ejTpw8xMTEAzJgxgw4dOtCmTRuGDBlCamoqAOfOnePOO++kTZs2tGnThk2bjBt25U7j3bJlS3r37m2fY+no0aP07duXdu3aceutt/L7778DcPz4cTp37kyHDh2YNGmSF65oKWmty81Pu3bttKeduZCie7/yg16x80+PH7u8iUuJ03EpcWaHISxi//79Ts+7zuqqZ+2cpbXWOiMrQ3ed1VV/vutzrbXWKRkpuuusrnr+nvlaa60T0hJ011ld9aL9i7TWxnuz66yu+vvfv9daax2TFFPsuZs1a6bj4+N1SkqK7tSpk3788ce11lqHhYU5bVu9enWttdarVq3Sbdu21QMHDtQJCQm6d+/e+sKFC0WeB9ALFiwwXlNGhu7cubOOjY3VWms9f/58/cADD2ittY6Pj7fv869//Uu///77Wmuthw0bpt99912ttdZZWVk6ISFBHz9+XPv7++udO3dqrbW+++679eefG9fp9ttv14cOHdJaa71lyxbdvXt3rbXWgwYN0nPmzNFaaz1t2jRdpUqVQq+Li9ewXXvpM9nyjdRnLqYA+EQJIjwk3OwQhPCI6667jmeffZZevXoRGhpKmzZtCAgo+uOqV69e9OrVC4A5c+bQv39/Dh48yJQpU6hRowZTp04lJMT5c8Df358hQ4YAcPDgQfbu3Ws/RnZ2NvXr1weMKq+JEyeSkJBAcnIyffr0AWDt2rV89tln9mOFhYVx8eJFl9N4Jycns2nTJu6++277+S9fNu50uXHjRhYtWgQYU5Y/++yzpb94HmT5BBFz0SgK+kKC+ObANwDcdd1dJkcirMhxGvlA/0Cn5yGBIU7Pw4LCnJ6Hh4Q7PXenM8VDDz3EQw89BMCECRNo1KgRAHXr1iUmJob69esTExNDnTp1nPZLTU1lzpw5rFy5kt69e7N48WK++OIL5s2bx5gxY5y2DQoKwt/fHzBqU1q2bMnmzZsLxDJ69Gi+++472rRpw+zZs1lXzNTlrqbxzsnJoXr16vYbBuVX2qnGvcnybRBnL6YS6O9n+TEQAO//+j7v//q+2WEI4RGxsbEA/Pnnn3zzzTfce++9gHFznjlz5gBGSWHw4MFO+7311ls8+eSTBAYGkpaWhlIKPz8/e7tBYVq0aEFcXJw9QWRmZrJv3z4AkpKSqF+/PpmZmcybN8++T48ePfjwww8Bo8Rx6dKlQo9frVo1mjRpwldffQUYCWnXrl2AMQOs45Tl5YXlE8SZi6nUqx6MXznMzp62ePhiFg9fbHYYQnjEkCFDuP766xk0aBAffPABNWrUAOC5555j9erVNGvWjNWrV/Pcc8/Z9zlz5gzbt2+3J41nnnmGTp06MWfOHPtd5wpTqVIlvv76a5599lnatGlDVFSUvdF58uTJ3HTTTfTq1Ytrr73Wvs/UqVP56aefuOGGG2jXrp09oRRm3rx5fPrpp7Rp04aWLVuyePFi+3E++OADOnToQGJiYskvlpdYfrrv//v4Z+qEBfPK8A4ePa4QVmf2dN+ioLKe7tvSJYjsnBxOnU8hIrx8DDrxtgV7F7Bgb8EbmAghRGlYupH67MU0MrNziKjtGwniw+1GXeg9re4pZkshhCiepRPEH/HGfWMjwquaHEnZWPbXZWaHIISwEEsniD/jkgF8poopJND6XXlF2dJal8vul77IjPZiS7dBnIhLona1IEIqWzoP2s3dPZe5u+eaHYawiKCgIM6fP2/KB5NwprXm/PnzBAWVbXd9S39yHjqTSLP6YWaHUWY++e0TAEa0HmFyJMIKGjVqxKlTp4iLizM7FIGRsHMHC5YVryYIpVRfYCrgD3yitX7Tm+dzlJSWyekLKfRqU7YX1Eyr719tdgjCQgIDA2nSpInZYQgTea2KSSnlD3wA9AOuB+5VSl3vrfPldygmAYDmDXynBBHoH0igf6DZYQghLMKbbRAdgSNa62Na6wxgPjC4mH08ZveJ8/gpxbUNqpfVKU03O3o2s6Nnmx2GEMIivFnF1BA46fD8FHBTUTucjE/mjW92Ur9GCNVCKhHor6gU4E/9GiFEhIcSFlLJ7R4Vvx2L59qG1akS5DvfqHOTw+io0abGIYSwBm8mCFef5AW6QyilxgJjbU8vTxjSdq8ng3jvQU8erUyFA/Gl2VE9YKluiaW+DhYk1yKPXIs8Lbx1YG8miFNAY4fnjYAz+TfSWk8HpgMopbZ7a06RikauhUGuQx65FnnkWuRRSnl2AjsH3myD2AY0U0o1UUpVAoYD33vxfEIIITzIayUIrXWWUupxYCVGN9eZWuui58IVQghRbnh1HITWehlQkgmCpnsrlgpIroVBrkMeuRZ55Frk8dq1KFf3gxBCCFF+WHouJiGEEKVXLhKEUqqvUuqgUuqIUuq54veomJRSJ5RSe5RS0bk9D5RSNZVSq5VSh22/a9iWK6XU+7Zrslsp1dbhOKNs2x9WSo0y6/WUhFJqplIqVim112GZx167Uqqd7doese1bbvv6FnItXlJKnba9N6KVUv0d1j1ve10HlVJ9HJa7/L+xdQz51XaNFtg6iZQ7SqnGSqmflFIHlFL7lFJP2pb73PuiiGth7vtCa23qD0YD9lGgKVAJ2AVcb3ZcXnqtJ4DwfMveAp6zPX4O+LftcX9gOcZ4kk7Ar7blNYFjtt81bI9rmP3a3HjttwFtgb3eeO3AVqCzbZ/lQD+zX3MJr8VLwHgX215v+5+oDDSx/a/4F/V/AywEhtsefwQ8YvZrLuQ61Afa2h5XBQ7ZXq/PvS+KuBamvi/KQwnC1Ck5yoHBwBzb4znAXxyWf6YNW4DqSqn6QB9gtdb6gtb6IrAa6FvWQZeU1vpn4EK+xR557bZ11bTWm7Xx7v/M4VjlTiHXojCDgfla68ta6+PAEYz/GZf/N7ZvyLcDX9v2d7yu5YrWOkZr/ZvtcRJwAGMGBp97XxRxLQpTJu+L8pAgXE3JUdSFqcg0sEoptUMZI8gB6mqtY8B4kwB1bMsLuy5Wul6eeu0NbY/zL69oHrdVnczMrVah5NeiFpCgtc7Kt7xcU0pFAjcCv+Lj74t81wJMfF+UhwTh1pQcFnGL1rotxgy3jymlbiti28Kuiy9cr5K+ditckw+Bq4EoIAZ4x7bc8tdCKRUKLALGaa0vFbWpi2VWvxamvi/KQ4Jwa0oOK9Ban7H9jgW+xSgOnrMVhbH9jrVtXth1sdL18tRrP2V7nH95haG1Pqe1ztZa5wAzMN4bUPJrEY9R9RKQb3m5pJQKxPhAnKe1/sa22CffF66uhdnvi/KQIHxiSg6lVBWlVNXcx0BvYC/Ga83tdTEKWGx7/D0w0tZzoxOQaCturwR6K6Vq2IqbvW3LKiKPvHbbuiSlVCdbXetIh2NVCLkfiDZ3Yrw3wLgWw5VSlZVSTYBmGA2vLv9vbHXtPwFDbfs7Xtdyxfa3+hQ4oLX+j8Mqn3tfFHYtTH9fmN16r/N6JxzCaH3/l9nxeOk1NsXoUbAL2Jf7OjHqBn8EDtt+17QtVxg3XDoK7AHaOxzrQYxGqSPAA2a/Njdf/5cYReRMjG85D3nytQPtbf88R4Fp2AaBlsefQq7F57bXutv2z1/fYft/2V7XQRx64RT2f2N7r221XaOvgMpmv+ZCrkMXjGqO3UC07ae/L74virgWpr4vZCS1EEIIl8pDFZMQQohySBKEEEIIlyRBCCGEcEkShBBCCJckQQghhHBJEoSokJRS65RSXr8nsVLqCdsMm/PyLW+vlHrf9ribUupmD54zUil1n6tzCVGWvHpHOSHKI6VUgM6bk6Y4j2L0MT/uuFBrvR3IvVl8NyAZ2OShGCKB+4AvXJxLiDIjJQjhNbZvwgeUUjNsc9yvUkoF29bZSwBKqXCl1Anb49FKqe+UUkuUUseVUo8rpZ5WSu1USm1RStV0OMUIpdQmpdRepVRH2/5VbJOabbPtM9jhuF8ppZYAq1zE+rTtOHuVUuNsyz7CGFz0vVLqqXzbd1NK/WCbWO3/gKeUMV//rUqp2kqpRbYYtimlbrHt85JSarpSahXwme36/KKU+s32k1sKeRO41Xa8p3LPZTtGTdv12W27Hq0djj3Tdl2PKaWecLgeS5VSu2yv7Z4r+6sKn2L2CEL5se4PxjfhLCDK9nwhMML2eB22kbBAOHDC9ng0xkjPqkBtIBH4P9u6dzEmMcvdf4bt8W3Y7q0AvO5wjuoYI0qr2I57Ctuo3HxxtsMYrVoFCMUY6X6jbd0J8t3Dw7a8G/CD7fFLOMzZj/HNv4vtcQTG9Am52+0Agm3PQ4Ag2+NmwPb8x3Zxrv8CL9oe3w5EOxx7E8b9AcKB80AgMCT3Otm2CzP7fSE/FedHqpiEtx3XWkfbHu/ASBrF+Ukbc+InKaUSgSW25XuA1g7bfQnG/RWUUtWUUtUx5uG5Qyk13rZNEMaHNNjuGeDifF2Ab7XWKQBKqW+AW4Gd7rxAF3oC16u8m5dVU7Z5uDDmxUmzPQ4EpimlooBsoLkbx+6C8aGP1nqtUqqWUirMtm6p1voycFkpFQvUxbhmU5RS/8ZIMr+U8jUJHyQJQnjbZYfH2UCw7XEWeVWcQUXsk+PwPAfn92z+eWJypzUeorU+6LhCKXUTkFJIjJ6+DaUf0NkhEeTGQL4YngLOAW1s+6S7ceyipm3Of60DtNaHlFLtMObneUMptUpr/Ypbr0L4PGmDEGY5gVG1A3kzTJbUPQBKqS4YM3smYszs+Xfb7JgopW504zg/A39RSoUoY6bdO4GSfNNOwqgSy7UKeDz3ia2E4EoYEKONqZzvx7hdpKvj5Y/1r7bjdgPidRH3UFBKNQBStdZzgSkYtzoVwi2SIIRZpgCPKKU2YdSZl8ZF2/4fYcyICjAZo+pmt1Jqr+15kbRxq8fZGDNd/gp8orUuSfXSEuDO3EZq4Amgva0heT9GI7Yr/wNGKaW2YFQv5ZYudgNZtoblp/Lt81LusTEas0dRtBuArUqpaIzZP18twesSPk5mcxVCCOGSlCCEEEK4JAlCCCGES5IghBBCuCQJQgghhEuSIIQQQrgkCUIIIYRLkiCEEEK4JAlCCCGES/8PMVbfKZDjDnEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.axhline(mi,label='ground truth',linestyle='--',color='red')\n",
    "for i in range(rep):\n",
    "    plt.plot(mi_list[i,:],color='steelblue')\n",
    "    for t in range(mi_list[i].shape[0]):\n",
    "        if (mi_list[0,t]>.9*mi):\n",
    "            plt.axvline(t,label='90% reached',linestyle=':',color='green')\n",
    "            break\n",
    "plt.xlim((0,mi_list[0].shape[0]))\n",
    "plt.ylim((0,mi*1.1))\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"MI estimate\")\n",
    "plt.legend()\n",
    "plt.savefig(fig_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
